| 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 // G4DNADoubleIonisationModel.cc
27 //
28 // Created at 2024/04/03 (Thu.)
29 // Author: Shogo OKADA @KEK-CRC (shogo.okada@kek.jp)
30 //
31 // Reference: J.Meesungnoen et. al, DOI: 10.1021/jp058037z
32 //
33
34 #include "G4DNADoubleIonisationModel.hh"
35
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UAtomicDeexcitation.hh"
39 #include "G4LossTableManager.hh"
40 #include "G4DNAChemistryManager.hh"
41 #include "G4DNAMolecularMaterial.hh"
42
43 #include "G4IonTable.hh"
44 #include "G4GenericIon.hh"
45 #include "G4DNARuddAngle.hh"
46 #include "G4DeltaAngle.hh"
47 #include "G4Exp.hh"
48
49 #include <sstream>
50
51 namespace {
52
53 G4DNAWaterIonisationStructure water_structure;
54
55 // parameters for rejection function
56 struct FuncParams {
57 G4double Bj_energy;
58 G4double alpha_const;
59 G4double beta_squared;
60 G4double velocity;
61 G4double correction_factor;
62 G4double wc;
63 G4double F1;
64 G4double F2;
65 G4double c;
66 };
67
68 //------------------------------------------------------------------------------
69 void setup_rejection_function(G4ParticleDefinition* pdef, const G4double ekin,
70 const G4int shell, FuncParams& par)
71 {
72
73 // Following values provided by M. Dingfelder (priv. comm)
74 const G4double Bj[5]
75 = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540.0 * eV };
76
77 // Data For Liquid Water from Dingfelder (Protons in Water)
78 G4double A1{1.02}, B1{82.0}, C1{0.45}, D1{-0.80}, E1{0.38}, A2{1.07},
79 B2{11.6}, // Value provided by M. Dingfelder (priv. comm)
80 C2{0.60}, D2{0.04}, alpha_const{0.64};
81
82 auto Bj_energy = Bj[shell];
83
84 if (shell == 4) {
85 alpha_const = 0.66;
86 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
87 A1 = 1.25; B1 = 0.5; C1 = 1.00; D1 = 1.00; E1 = 3.00;
88 A2 = 1.10; B2 = 1.30; C2 = 1.00; D2 = 0.00;
89 // The following cases are provided by M. Dingfelder (priv. comm)
90 Bj_energy = water_structure.IonisationEnergy(shell);
91 }
92
93 const auto mass = pdef->GetPDGMass();
94 const auto tau = ekin * electron_mass_c2 / mass;
95 const auto A_ion = pdef->GetAtomicMass();
96
97 G4double v2;
98 G4double beta2;
99
100 constexpr G4double Ry = 13.6 * eV;
101 constexpr G4double xxx = 5.447761194E-02 * MeV;
102
103 if (tau < xxx) {
104 v2 = tau / Bj_energy;
105 beta2 = 2.0 * tau / electron_mass_c2;
106 } else {
107 // Relativistic
108 v2 = (0.5 * electron_mass_c2 / Bj_energy)
109 * (1.0 - (1.0 / std::pow((1.0 + (tau / electron_mass_c2)), 2.0)));
110 beta2 = 1.0 - 1.0 / std::pow((1.0 + (tau / electron_mass_c2 / A_ion)), 2.0);
111 }
112
113 const auto v = std::sqrt(v2);
114 const auto wc = 4.0 * v2 - 2.0 * v - (Ry / (4.0 * Bj_energy));
115
116 const auto L1 = (C1 * std::pow(v, D1)) / (1.0 + E1 * std::pow(v, (D1 + 4.0)));
117 const auto L2 = C2 * std::pow(v, D2);
118 const auto H1 = (A1 * G4Log(1.0 + v2)) / (v2 + (B1 / v2));
119 const auto H2 = (A2 / v2) + (B2 /(v2 * v2));
120 const auto F1 = L1 + H1;
121 const auto F2 = (L2 * H2) / (L2 + H2);
122
123 // ZF. generalized & relativistic version
124 G4double max_energy;
125
126 if (ekin <= 0.1 * mass) {
127 // maximum kinetic energy , non relativistic
128 max_energy = 4.0 * (electron_mass_c2 / mass) * ekin;
129 } else {
130 // relativistic
131 auto gamma = 1.0 / std::sqrt(1.0 - beta2);
132 max_energy = 2.0 * electron_mass_c2 * (gamma * gamma - 1.0)
133 / (1.0 + 2.0 * gamma * (electron_mass_c2 / mass)
134 + std::pow(electron_mass_c2 / mass, 2.0));
135 }
136
137 const auto wmax = max_energy / Bj_energy;
138 auto c = wmax * (F2 * wmax+ F1 * (2.0 + wmax))
139 / (2.0 * (1.0 + wmax) * (1.0 + wmax));
140 c = 1.0 / c; // manual calculus leads to c = 1 / c
141
142 par.Bj_energy = Bj_energy;
143 par.alpha_const = alpha_const;
144 par.beta_squared = beta2;
145 par.velocity = v;
146 par.correction_factor = 1.0;
147 par.wc = wc;
148 par.F1 = F1;
149 par.F2 = F2;
150 par.c = c;
151
152 }
153
154 //------------------------------------------------------------------------------
155 G4double rejection_function(G4ParticleDefinition* pdef, const G4int shell,
156 const FuncParams& par, G4double proposed_ws)
157 {
158
159 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1.0 };
160
161 proposed_ws /= par.Bj_energy;
162
163 auto rejection_term = 1.0 + G4Exp(par.alpha_const * (proposed_ws - par.wc)
164 / par.velocity);
165 rejection_term = (1.0 / rejection_term) * par.correction_factor * Gj[shell];
166
167 if (pdef == G4Proton::ProtonDefinition()) {
168
169 // for protons
170 return rejection_term;
171
172 } else if (pdef->GetAtomicMass() > 4) {
173
174 // for carbon ions
175 auto Z = pdef->GetAtomicNumber();
176 auto x = 100.0 * std::sqrt(par.beta_squared) / std::pow(Z, 0.6666667);
177 auto zeff = Z * (1.0 - G4Exp(x * (-1.316 + x * (0.112 - 0.0650 * x))));
178
179 rejection_term *= (zeff * zeff);
180
181 return rejection_term;
182
183 }
184
185 // for alpha particles
186 auto zeff = pdef->GetPDGCharge() / eplus + pdef->GetLeptonNumber();
187 rejection_term *= (zeff * zeff);
188
189 return rejection_term;
190
191 }
192
193 //------------------------------------------------------------------------------
194 G4double proposed_sampled_energy(const FuncParams& par)
195 {
196 const auto rval = G4UniformRand();
197
198 auto proposed_ws = par.c * (par.F1 * par.F1 * par.c
199 + 2.0 * rval * (par.F2 - par.F1));
200
201 proposed_ws = -par.F1 * par.c + 2.0 * rval + std::sqrt(proposed_ws);
202
203 proposed_ws /= (par.c * (par.F1 + par.F2) - 2.0 * rval);
204
205 proposed_ws *= par.Bj_energy;
206
207 return proposed_ws;
208 }
209
210 } // end of anonymous namespace
211
212 //==============================================================================
213
214 // constructor
215 G4DNADoubleIonisationModel::G4DNADoubleIonisationModel(
216 const G4ParticleDefinition*, const G4String& model_name)
217 : G4VEmModel(model_name),
218 is_initialized_(false)
219 {
220 water_density_ = nullptr;
221
222 model_elow_tab_[1] = 100 * eV;
223 model_elow_tab_[4] = 1.0 * keV;
224 model_elow_tab_[5] = 0.5 * MeV; // For A = 3 or above, limit is MeV/uma
225
226 verbose_level_ = 0;
227
228 // Define default angular generator
229 SetAngularDistribution(new G4DNARuddAngle());
230
231 // Mark this model as "applicable" for atomic deexcitation
232 SetDeexcitationFlag(true);
233 atom_deex_ = nullptr;
234 particle_change_ = nullptr;
235
236 // Selection of stationary mode
237 stat_code_ = false;
238
239 // True if use champion alpha parameter
240 use_champion_param_ = false;
241
242 // Double-ionization energy
243 energy_threshold_ = 40.0 * eV;
244 }
245
246 //------------------------------------------------------------------------------
247 G4DNADoubleIonisationModel::~G4DNADoubleIonisationModel()
248 {
249 for (const auto& x : xs_tab_) {
250 G4DNACrossSectionDataSet* table = x.second;
251 if (table) { delete table; }
252 }
253 }
254
255 //------------------------------------------------------------------------------
256 void G4DNADoubleIonisationModel::Initialise(
257 const G4ParticleDefinition* particle, const G4DataVector&)
258 {
259 if (verbose_level_ > 3) {
260 G4cout << "Calling G4DNADoubleIonisationModel::Initialise()" << G4endl;
261 }
262
263 proton_def_ = G4Proton::ProtonDefinition();
264 alpha_def_ = G4DNAGenericIonsManager::Instance()->GetIon("alpha++");
265 carbon_def_ = G4IonTable::GetIonTable()->GetIon(6, 12);
266
267 constexpr G4double kScaleFactor = 1.0 * m * m;
268
269 mioni_manager_ = new G4DNAMultipleIonisationManager();
270
271 G4double Z{0.0}, A{0.0};
272 G4String alpha_param_file{"dna/multipleionisation_alphaparam_champion.dat"};
273
274 if (particle == proton_def_) {
275
276 // *************************************************************************
277 // for protons
278 const auto& proton = proton_def_->GetParticleName();
279 elow_tab_[proton] = model_elow_tab_[1];
280 eupp_tab_[proton] = 3.0 * MeV;
281
282 // load cross-section data for single ionization process
283 auto xs_proton = new G4DNACrossSectionDataSet(
284 new G4LogLogInterpolation, eV, kScaleFactor);
285 xs_proton->LoadData("dna/sigma_ionisation_p_rudd");
286 xs_tab_[proton] = xs_proton;
287
288 // set energy limits
289 SetLowEnergyLimit(elow_tab_[proton]);
290 SetHighEnergyLimit(eupp_tab_[proton]);
291
292 if (!use_champion_param_) {
293 alpha_param_file = "dna/multipleionisation_alphaparam_p.dat";
294 }
295
296 Z = static_cast<G4double>(proton_def_->GetAtomicNumber());
297 A = static_cast<G4double>(proton_def_->GetAtomicMass());
298
299 } else if (particle == alpha_def_) {
300
301 //**************************************************************************
302 // for alpha particles
303 const auto& alpha = alpha_def_->GetParticleName();
304 elow_tab_[alpha] = model_elow_tab_[4];
305 eupp_tab_[alpha] = 23.0 * MeV;
306
307 // load cross-section data for single ionization process
308 auto xs_alpha = new G4DNACrossSectionDataSet(
309 new G4LogLogInterpolation, eV, kScaleFactor);
310 xs_alpha->LoadData("dna/sigma_ionisation_alphaplusplus_rudd");
311 xs_tab_[alpha] = xs_alpha;
312
313 // set energy limits
314 SetLowEnergyLimit(elow_tab_[alpha]);
315 SetHighEnergyLimit(eupp_tab_[alpha]);
316
317 if (!use_champion_param_) {
318 alpha_param_file = "dna/multipleionisation_alphaparam_alphaplusplus.dat";
319 }
320
321 Z = static_cast<G4double>(alpha_def_->GetAtomicNumber());
322 A = static_cast<G4double>(alpha_def_->GetAtomicMass());
323
324 } else if (particle == G4GenericIon::GenericIonDefinition()) {
325
326 // *************************************************************************
327 // for carbon ions
328 const auto& carbon = carbon_def_->GetParticleName();
329 elow_tab_[carbon] = model_elow_tab_[5] * carbon_def_->GetAtomicMass();
330 eupp_tab_[carbon] = 120.0 * MeV;
331
332 // load cross-section data for single ionization process
333 auto xs_carbon = new G4DNACrossSectionDataSet(
334 new G4LogLogInterpolation, eV, kScaleFactor);
335 xs_carbon->LoadData("dna/sigma_ionisation_c_rudd");
336 xs_tab_[carbon] = xs_carbon;
337
338 // set energy limits
339 SetLowEnergyLimit(elow_tab_[carbon]);
340 SetHighEnergyLimit(eupp_tab_[carbon]);
341
342 if (!use_champion_param_) {
343 alpha_param_file = "dna/multipleionisation_alphaparam_c.dat";
344 }
345
346 Z = static_cast<G4double>(carbon_def_->GetAtomicNumber());
347 A = static_cast<G4double>(carbon_def_->GetAtomicMass());
348
349 }
350
351 // load alpha parameter
352 mioni_manager_->LoadAlphaParam(alpha_param_file, Z, A);
353
354 if (verbose_level_ > 0) {
355 G4cout << "G4DNADoubleIonisationModel is initialized " << G4endl
356 << "Energy range: "
357 << LowEnergyLimit() / eV << " eV - "
358 << HighEnergyLimit() / keV << " keV for "
359 << particle->GetParticleName()
360 << G4endl;
361 }
362
363 water_density_ = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(
364 G4Material::GetMaterial("G4_WATER"));
365
366 atom_deex_ = G4LossTableManager::Instance()->AtomDeexcitation();
367
368 if (is_initialized_) { return; }
369
370 particle_change_ = GetParticleChangeForGamma();
371 is_initialized_ = true;
372 }
373
374 //------------------------------------------------------------------------------
375 G4double G4DNADoubleIonisationModel::GetLowEnergyLimit(const G4String& pname)
376 {
377 G4double elim{0.0};
378
379 EnergyLimitTable::iterator itr = elow_tab_.find(pname);
380 if (itr != elow_tab_.end()) { elim = itr->second; }
381
382 return elim;
383 }
384
385 //------------------------------------------------------------------------------
386 G4double G4DNADoubleIonisationModel::GetUppEnergyLimit(const G4String& pname)
387 {
388 G4double elim{0.0};
389
390 EnergyLimitTable::iterator itr = eupp_tab_.find(pname);
391 if (itr != eupp_tab_.end()) { elim = itr->second; }
392
393 return elim;
394 }
395
396 //------------------------------------------------------------------------------
397 G4double G4DNADoubleIonisationModel::CrossSectionPerVolume(
398 const G4Material* material, const G4ParticleDefinition* pdef,
399 G4double ekin, G4double, G4double)
400 {
401
402 if (verbose_level_ > 3) {
403 G4cout << "Calling G4DNADoubleIonisationModel::CrossSectionPerVolume()"
404 << G4endl;
405 }
406
407 // Calculate total cross section for model
408
409 if (pdef != proton_def_ && pdef != alpha_def_ && pdef != carbon_def_) {
410 return 0.0;
411 }
412
413 static G4double water_dens = (*water_density_)[material->GetIndex()];
414
415 const auto& pname = pdef->GetParticleName();
416
417 const auto low_energy_lim = GetLowEnergyLimit(pname);
418 const auto upp_energy_lim = GetUppEnergyLimit(pname);
419
420 G4double sigma{0.0};
421 if (ekin <= upp_energy_lim) {
422
423 if (ekin < low_energy_lim) { ekin = low_energy_lim; }
424
425 CrossSectionDataTable::iterator pos = xs_tab_.find(pname);
426 if (pos == xs_tab_.end()) {
427 G4Exception("G4DNADoubleIonisationModel::CrossSectionPerVolume",
428 "em0002", FatalException,
429 "Model not applicable to particle type.");
430 }
431
432 G4DNACrossSectionDataSet* table = pos->second;
433 if (table != nullptr) {
434 const auto a = mioni_manager_->GetAlphaParam(ekin);
435 sigma = table->FindValue(ekin) * a;
436 }
437
438 }
439
440 if (verbose_level_ > 2) {
441
442 std::stringstream msg;
443
444 msg << "----------------------------------------------------------------\n";
445 msg << " G4DNADoubleIonisationModel - XS INFO START\n";
446 msg << " - Kinetic energy(eV): " << ekin/eV << ", Particle : "
447 << pdef->GetParticleName() << "\n";
448 msg << " - Cross section per water molecule (cm^2): "
449 << sigma / cm / cm << "\n";
450 msg << " - Cross section per water molecule (cm^-1): "
451 << sigma * water_dens / (1.0 / cm) << "\n";
452 msg << " G4DNADoubleIonisationModel - XS INFO END\n";
453 msg << "----------------------------------------------------------------\n";
454
455 G4cout << msg.str() << G4endl;
456
457 }
458
459 return (sigma * water_dens);
460
461 }
462
463 //------------------------------------------------------------------------------
464 G4double G4DNADoubleIonisationModel::GenerateSecondaries(
465 std::vector<G4DynamicParticle*>* vsec, const G4MaterialCutsCouple* couple,
466 const G4DynamicParticle* particle, G4int ioni_shell,
467 G4double& theta, G4double& phi, G4double& shell_energy)
468 {
469 auto pdef = particle->GetDefinition();
470
471 // get kinetic energy for a parent particle
472 auto ekin1 = particle->GetKineticEnergy();
473
474 // sample kinetic energy for a secondary electron
475 auto ekin2 = RandomizeEjectedElectronEnergy(pdef, ekin1, ioni_shell);
476
477 // sample momentum direction for a secondary electron
478 auto sample_electron_direction = [this](
479 const G4DynamicParticle* dp, G4double _ekin2, G4int _Z, G4int _ioni_shell,
480 const G4MaterialCutsCouple* mcc, G4double& _theta, G4double& _phi) {
481
482 G4ThreeVector locdir;
483
484 if (_theta > 0.0) {
485
486 auto costh = std::cos(_theta);
487 auto sinth = std::sqrt((1.0 - costh) * (1.0 + costh));
488 locdir.set(sinth * std::cos(_phi), sinth * std::sin(_phi), costh);
489 locdir.rotateUz(dp->GetMomentumDirection());
490
491 } else {
492
493 locdir = GetAngularDistribution()->SampleDirectionForShell(
494 dp, _ekin2, _Z, _ioni_shell, mcc->GetMaterial());
495 _theta = locdir.theta();
496 _phi = locdir.phi();
497
498 }
499
500 return locdir;
501 };
502
503 constexpr G4int Z = 8;
504 auto delta_dir = sample_electron_direction(
505 particle, ekin2, Z, ioni_shell, couple, theta, phi);
506
507 // generate a secondary electron and put it into the stack
508 auto dp = new G4DynamicParticle(G4Electron::Electron(), delta_dir, ekin2);
509 vsec->push_back(dp);
510
511 if (!atom_deex_ || ioni_shell != 4) { return ekin2; }
512
513 // ***************************************************************************
514 // Only atomic deexcitation from K shell is considered
515
516 constexpr auto k_shell = G4AtomicShellEnumerator(0);
517 const auto shell = atom_deex_->GetAtomicShell(Z, k_shell);
518
519 // get number of secondary electrons in the stack
520 // before processing atomic deescitation
521 const auto num_sec_init = vsec->size();
522
523 // perform atomic deexcitation process
524 atom_deex_->GenerateParticles(vsec, shell, Z, 0, 0);
525
526 // get number of secondary electrons in the stack
527 // after processing atomic deescitation
528 const auto num_sec_final = vsec->size();
529
530 if (num_sec_final == num_sec_init) { return ekin2; }
531
532 for (auto i = num_sec_init; i < num_sec_final; i++) {
533
534 auto e = ((*vsec)[i])->GetKineticEnergy();
535
536 // Check if there is enough residual energy
537 if (shell_energy < e) {
538
539 // Invalid secondary: not enough energy to create it!
540 // Keep its energy in the local deposit
541 delete (*vsec)[i];
542 (*vsec)[i] = 0;
543
544 continue;
545
546 }
547
548 // Ok, this is a valid secondary: keep it
549 shell_energy -= e;
550 }
551
552 // ***************************************************************************
553
554 return ekin2;
555 }
556
557 //------------------------------------------------------------------------------
558 void G4DNADoubleIonisationModel::SampleSecondaries(
559 std::vector<G4DynamicParticle*>* vsec, const G4MaterialCutsCouple* couple,
560 const G4DynamicParticle* particle, G4double, G4double)
561 {
562
563 if (verbose_level_ > 3) {
564 G4cout << "Calling SampleSecondaries() of G4DNADoubleIonisationModel"
565 << G4endl;
566 }
567
568 // get the definition for this parent particle
569 auto pdef = particle->GetDefinition();
570
571 // get kinetic energy
572 auto ekin = particle->GetKineticEnergy();
573
574 // get particle name
575 const auto& pname = pdef->GetParticleName();
576
577 // get energy limits
578 const auto low_energy_lim = GetLowEnergyLimit(pname);
579
580 // ***************************************************************************
581 // stop the transportation process of this parent particle
582 // if its kinetic energy is below the lower limit
583 if (ekin < low_energy_lim) {
584 particle_change_->SetProposedKineticEnergy(0.0);
585 particle_change_->ProposeTrackStatus(fStopAndKill);
586 particle_change_->ProposeLocalEnergyDeposit(ekin);
587 return;
588 }
589 // ***************************************************************************
590
591 constexpr G4int kNumSecondaries = 2;
592 constexpr G4double kDeltaTheta = pi;
593
594 G4int ioni_shell[kNumSecondaries];
595 G4double shell_energy[kNumSecondaries];
596
597 const auto scale_param = mioni_manager_->GetAlphaParam(ekin);
598 G4double tot_ioni_energy{0.0};
599 for (G4int i = 0; i < kNumSecondaries; i++) {
600 ioni_shell[i] = RandomSelect(ekin, scale_param, pname);
601 shell_energy[i] = ::water_structure.IonisationEnergy(ioni_shell[i]);
602 tot_ioni_energy += shell_energy[i];
603 }
604
605 if (ekin < tot_ioni_energy || tot_ioni_energy < energy_threshold_) {
606 return;
607 }
608
609 // generate secondary electrons
610 G4double theta{0.0}, phi{0.0}, tot_ekin2{0.0};
611 for (G4int i = 0; i < kNumSecondaries; i++) {
612 tot_ekin2 += GenerateSecondaries(vsec, couple, particle, ioni_shell[i],
613 theta, phi, shell_energy[i]);
614 theta += kDeltaTheta;
615 }
616
617 // This should never happen
618 if (mioni_manager_->CheckShellEnergy(eDoubleIonisedMolecule, shell_energy)) {
619 G4Exception("G4DNADoubleIonisatioModel::SampleSecondaries()",
620 "em2050", FatalException, "Negative local energy deposit");
621 }
622
623 // ***************************************************************************
624 // update kinematics for this parent particle
625 const auto primary_dir = particle->GetMomentumDirection();
626 particle_change_->ProposeMomentumDirection(primary_dir);
627
628 const auto scattered_energy = ekin - tot_ioni_energy - tot_ekin2;
629
630 // update total amount of shell energy
631 tot_ioni_energy = shell_energy[0] + shell_energy[1];
632
633 if (stat_code_) {
634 particle_change_->SetProposedKineticEnergy(ekin);
635 particle_change_->ProposeLocalEnergyDeposit(ekin - scattered_energy);
636 } else {
637 particle_change_->SetProposedKineticEnergy(scattered_energy);
638 particle_change_->ProposeLocalEnergyDeposit(tot_ioni_energy);
639 }
640
641 // ***************************************************************************
642 // generate double-ionized water molecules (H2O^2+)
643 const auto the_track = particle_change_->GetCurrentTrack();
644 mioni_manager_->CreateMultipleIonisedWaterMolecule(
645 eDoubleIonisedMolecule, ioni_shell, the_track);
646 // ***************************************************************************
647
648 }
649
650 //------------------------------------------------------------------------------
651 G4double G4DNADoubleIonisationModel::RandomizeEjectedElectronEnergy(
652 G4ParticleDefinition* pdef, G4double ekin, G4int shell)
653 {
654
655 //
656 // based on RandomizeEjectedElectronEnergy()
657 // of G4DNARuddIonisationExtendedModel
658 //
659
660 ::FuncParams par;
661 ::setup_rejection_function(pdef, ekin, shell, par);
662
663 // calculate maximum value
664 G4double emax{0.0}, val;
665 for (G4double en = 0.0; en < 20.0; en += 1.0) {
666 val = ::rejection_function(pdef, shell, par, en);
667 if (val <= emax) { continue; }
668 emax = val;
669 }
670
671 G4double proposed_energy, rand;
672 do {
673 // Proposed energy by inverse function sampling
674 proposed_energy = ::proposed_sampled_energy(par);
675 rand = G4UniformRand() * emax;
676 val = ::rejection_function(pdef, shell, par, proposed_energy);
677 } while (rand > val);
678
679 return proposed_energy;
680 }
681
682 //------------------------------------------------------------------------------
683 G4int G4DNADoubleIonisationModel::RandomSelect(
684 G4double ekin, G4double scale_param, const G4String& pname)
685 {
686
687 //
688 // based on RandomSelect() of G4DNARuddIonisationExtendedModel
689 //
690
691 // Retrieve data table corresponding to the current particle type
692 CrossSectionDataTable::iterator pos = xs_tab_.find(pname);
693
694 if (pos == xs_tab_.end()) {
695 G4Exception("G4DNADoubleIonisationModel::RandomSelect", "em0002",
696 FatalException, "Model not applicable to particle type.");
697 }
698
699 G4DNACrossSectionDataSet* table = pos->second;
700
701 if (table != nullptr) {
702
703 // get total number of energy level
704 const auto num_component = table->NumberOfComponents();
705
706 auto* valuesBuffer = new G4double[num_component];
707
708 auto shell = num_component;
709 G4double value = 0.0;
710
711 while (shell > 0) {
712 shell--;
713 valuesBuffer[shell] = table->GetComponent((G4int)shell)->FindValue(ekin)
714 * scale_param;
715 value += valuesBuffer[shell];
716 }
717
718 value *= G4UniformRand();
719
720 shell = num_component;
721
722 while (shell > 0) {
723 shell--;
724 if (valuesBuffer[shell] > value) {
725 delete [] valuesBuffer;
726 return (G4int)shell;
727 }
728 value -= valuesBuffer[shell];
729 }
730
731 delete [] valuesBuffer;
732 }
733
734 return 0;
735 }
736