| 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 // G4DNATripleIonisationModel.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 "G4DNATripleIonisationModel.hh"
35
36 #include "G4PhysicalConstants.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4UAtomicDeexcitation.hh"
39 #include "G4LossTableManager.hh"
40
41 #include "G4SystemOfUnits.hh"
42 #include "G4DNAMolecularMaterial.hh"
43 #include "G4IonTable.hh"
44 #include "G4GenericIon.hh"
45 #include "G4DNARuddAngle.hh"
46
47 #include <sstream>
48
49 namespace {
50
51 G4DNAWaterIonisationStructure water_structure;
52
53 } // end of anonymous namespace
54
55 //==============================================================================
56
57 // constructor
58 G4DNATripleIonisationModel::G4DNATripleIonisationModel(
59 const G4ParticleDefinition* p, const G4String& model_name)
60 : G4DNADoubleIonisationModel(p, model_name)
61 {
62 // Triple-ionisation energy
63 energy_threshold_ = 65.0 * eV;
64 }
65
66 //------------------------------------------------------------------------------
67 void G4DNATripleIonisationModel::Initialise(
68 const G4ParticleDefinition* particle, const G4DataVector&)
69 {
70 if (verbose_level_ > 3) {
71 G4cout << "Calling G4DNATripleIonisationModel::Initialise()" << G4endl;
72 }
73
74 proton_def_ = G4Proton::ProtonDefinition();
75 alpha_def_ = G4DNAGenericIonsManager::Instance()->GetIon("alpha++");
76 carbon_def_ = G4IonTable::GetIonTable()->GetIon(6, 12);
77
78 constexpr G4double kScaleFactor = 1.0 * m * m;
79
80 mioni_manager_ = new G4DNAMultipleIonisationManager();
81
82 G4double Z{0.0}, A{0.0};
83 G4String alpha_param_file{"dna/multipleionisation_alphaparam_champion.dat"};
84
85 if (particle == proton_def_) {
86
87 // *************************************************************************
88 // for protons
89 const auto& proton = proton_def_->GetParticleName();
90 elow_tab_[proton] = model_elow_tab_[1];
91 eupp_tab_[proton] = 3.0 * MeV;
92
93 // load cross-section data for single ionization process
94 auto xs_proton = new G4DNACrossSectionDataSet(
95 new G4LogLogInterpolation, eV, kScaleFactor);
96 xs_proton->LoadData("dna/sigma_ionisation_p_rudd");
97 xs_tab_[proton] = xs_proton;
98
99 // set energy limits
100 SetLowEnergyLimit(elow_tab_[proton]);
101 SetHighEnergyLimit(eupp_tab_[proton]);
102
103 if (!use_champion_param_) {
104 alpha_param_file = "dna/multipleionisation_alphaparam_p.dat";
105 }
106
107 Z = static_cast<G4double>(proton_def_->GetAtomicNumber());
108 A = static_cast<G4double>(proton_def_->GetAtomicMass());
109
110 } else if (particle == alpha_def_) {
111
112 //**************************************************************************
113 // for alpha particles
114 const auto& alpha = alpha_def_->GetParticleName();
115 elow_tab_[alpha] = model_elow_tab_[4];
116 eupp_tab_[alpha] = 23.0 * MeV;
117
118 // load cross-section data for single ionization process
119 auto xs_alpha = new G4DNACrossSectionDataSet(
120 new G4LogLogInterpolation, eV, kScaleFactor);
121 xs_alpha->LoadData("dna/sigma_ionisation_alphaplusplus_rudd");
122 xs_tab_[alpha] = xs_alpha;
123
124 // set energy limits
125 SetLowEnergyLimit(elow_tab_[alpha]);
126 SetHighEnergyLimit(eupp_tab_[alpha]);
127
128 if (!use_champion_param_) {
129 alpha_param_file = "dna/multipleionisation_alphaparam_alphaplusplus.dat";
130 }
131
132 Z = static_cast<G4double>(alpha_def_->GetAtomicNumber());
133 A = static_cast<G4double>(alpha_def_->GetAtomicMass());
134
135 } else if (particle == G4GenericIon::GenericIonDefinition()) {
136
137 // *************************************************************************
138 // for carbon ions
139 const auto& carbon = carbon_def_->GetParticleName();
140 elow_tab_[carbon] = model_elow_tab_[5] * carbon_def_->GetAtomicMass();
141 eupp_tab_[carbon] = 120.0 * MeV;
142
143 // load cross-section data for single ionization process
144 auto xs_carbon = new G4DNACrossSectionDataSet(
145 new G4LogLogInterpolation, eV, kScaleFactor);
146 xs_carbon->LoadData("dna/sigma_ionisation_c_rudd");
147 xs_tab_[carbon] = xs_carbon;
148
149 // set energy limits
150 SetLowEnergyLimit(elow_tab_[carbon]);
151 SetHighEnergyLimit(eupp_tab_[carbon]);
152
153 if (!use_champion_param_) {
154 alpha_param_file = "dna/multipleionisation_alphaparam_c.dat";
155 }
156
157 Z = static_cast<G4double>(carbon_def_->GetAtomicNumber());
158 A = static_cast<G4double>(carbon_def_->GetAtomicMass());
159
160 }
161
162 // load alpha parameter
163 mioni_manager_->LoadAlphaParam(alpha_param_file, Z, A);
164
165 if (verbose_level_ > 0) {
166 G4cout << "G4DNATripleIonisationModel is initialized " << G4endl
167 << "Energy range: "
168 << LowEnergyLimit() / eV << " eV - "
169 << HighEnergyLimit() / keV << " keV for "
170 << particle->GetParticleName()
171 << G4endl;
172 }
173
174 water_density_ = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(
175 G4Material::GetMaterial("G4_WATER"));
176
177 atom_deex_ = G4LossTableManager::Instance()->AtomDeexcitation();
178
179 if (is_initialized_) { return; }
180
181 particle_change_ = GetParticleChangeForGamma();
182 is_initialized_ = true;
183 }
184
185 //------------------------------------------------------------------------------
186 G4double G4DNATripleIonisationModel::CrossSectionPerVolume(
187 const G4Material* material, const G4ParticleDefinition* pdef,
188 G4double ekin, G4double, G4double)
189 {
190
191 if (verbose_level_ > 3) {
192 G4cout << "Calling G4DNATripleIonisationModel::CrossSectionPerVolume()"
193 << G4endl;
194 }
195
196 // Calculate total cross section for model
197
198 if (pdef != proton_def_ && pdef != alpha_def_ && pdef != carbon_def_) {
199 return 0.0;
200 }
201
202 static G4double water_dens = (*water_density_)[material->GetIndex()];
203
204 const auto& pname = pdef->GetParticleName();
205
206 const auto low_energy_lim = GetLowEnergyLimit(pname);
207 const auto upp_energy_lim = GetUppEnergyLimit(pname);
208
209 G4double sigma{0.0};
210 if (ekin <= upp_energy_lim) {
211
212 if (ekin < low_energy_lim) { ekin = low_energy_lim; }
213
214 CrossSectionDataTable::iterator pos = xs_tab_.find(pname);
215 if (pos == xs_tab_.end()) {
216 G4Exception("G4DNATripleIonisationModel::CrossSectionPerVolume",
217 "em0002", FatalException,
218 "Model not applicable to particle type.");
219 }
220
221 G4DNACrossSectionDataSet* table = pos->second;
222 if (table != nullptr) {
223 const auto a = mioni_manager_->GetAlphaParam(ekin);
224 sigma = table->FindValue(ekin) * a * a;
225 }
226
227 }
228
229 if (verbose_level_ > 2) {
230
231 std::stringstream msg;
232
233 msg << "----------------------------------------------------------------\n";
234 msg << " G4DNATripleIonisationModel - XS INFO START\n";
235 msg << " - Kinetic energy(eV): " << ekin/eV << ", Particle : "
236 << pdef->GetParticleName() << "\n";
237 msg << " - Cross section per water molecule (cm^2): "
238 << sigma / cm / cm << "\n";
239 msg << " - Cross section per water molecule (cm^-1): "
240 << sigma * water_dens / (1.0 / cm) << "\n";
241 msg << " G4DNATripleIonisationModel - XS INFO END\n";
242 msg << "----------------------------------------------------------------\n";
243
244 G4cout << msg.str() << G4endl;
245
246 }
247
248 return (sigma * water_dens);
249
250 }
251
252 //------------------------------------------------------------------------------
253 void G4DNATripleIonisationModel::SampleSecondaries(
254 std::vector<G4DynamicParticle*>* vsec, const G4MaterialCutsCouple* couple,
255 const G4DynamicParticle* particle, G4double, G4double)
256 {
257
258 if (verbose_level_ > 3) {
259 G4cout << "Calling SampleSecondaries() of G4DNATripleIonisationModel"
260 << G4endl;
261 }
262
263 // get the definition for this parent particle
264 auto pdef = particle->GetDefinition();
265
266 // get kinetic energy
267 auto ekin = particle->GetKineticEnergy();
268
269 // get particle name
270 const auto& pname = pdef->GetParticleName();
271
272 // get energy limits
273 const auto low_energy_lim = GetLowEnergyLimit(pname);
274
275 // ***************************************************************************
276 // stop the transportation process of this parent particle
277 // if its kinetic energy is below the lower limit
278 if (ekin < low_energy_lim) {
279 particle_change_->SetProposedKineticEnergy(0.0);
280 particle_change_->ProposeTrackStatus(fStopAndKill);
281 particle_change_->ProposeLocalEnergyDeposit(ekin);
282 return;
283 }
284 // ***************************************************************************
285
286 constexpr G4int kNumSecondaries = 3;
287 constexpr G4double kDeltaTheta = pi * 0.666666667;
288
289 G4int ioni_shell[kNumSecondaries] = {0, 0, 0};
290 G4double shell_energy[kNumSecondaries];
291
292 auto scale_param = mioni_manager_->GetAlphaParam(ekin);
293 scale_param *= scale_param;
294
295 G4bool is_continue{true};
296 while (1) {
297 ioni_shell[0] = RandomSelect(ekin, scale_param, pname);
298 ioni_shell[1] = RandomSelect(ekin, scale_param, pname);
299 ioni_shell[2] = RandomSelect(ekin, scale_param, pname);
300 is_continue = (ioni_shell[0] == ioni_shell[1] &&
301 ioni_shell[1] == ioni_shell[2]);
302 if (!is_continue) { break; }
303 }
304
305 G4double tot_ioni_energy{0.0};
306 for (int i = 0; i < kNumSecondaries; i++) {
307 shell_energy[i] = ::water_structure.IonisationEnergy(ioni_shell[i]);
308 tot_ioni_energy += shell_energy[i];
309 }
310
311 if (ekin < tot_ioni_energy || tot_ioni_energy < energy_threshold_) {
312 return;
313 }
314
315 // generate secondary electrons
316 G4double theta{0.0}, phi{0.0}, tot_ekin2{0.0};
317 for (int i = 0; i < kNumSecondaries; i++) {
318 tot_ekin2 += GenerateSecondaries(vsec, couple, particle, ioni_shell[i],
319 theta, phi, shell_energy[i]);
320 theta += kDeltaTheta;
321 }
322
323 // This should never happen
324 if (mioni_manager_->CheckShellEnergy(eTripleIonisedMolecule, shell_energy)) {
325 G4Exception("G4DNATripleIonisatioModel::SampleSecondaries()",
326 "em2050", FatalException, "Negative local energy deposit");
327 }
328
329 // ***************************************************************************
330 // update kinematics for this parent particle
331 const auto primary_dir = particle->GetMomentumDirection();
332 particle_change_->ProposeMomentumDirection(primary_dir);
333
334 const auto scattered_energy = ekin - tot_ioni_energy - tot_ekin2;
335
336 // update total amount of shell energy
337 tot_ioni_energy = shell_energy[0] + shell_energy[1] + shell_energy[2];
338
339 if (stat_code_) {
340 particle_change_->SetProposedKineticEnergy(ekin);
341 particle_change_->ProposeLocalEnergyDeposit(
342 ekin - scattered_energy);
343 } else {
344 particle_change_->SetProposedKineticEnergy(scattered_energy);
345 particle_change_->ProposeLocalEnergyDeposit(tot_ioni_energy);
346 }
347
348 // ***************************************************************************
349 // generate triple-ionized water molecules (H2O^3+)
350 const auto the_track = particle_change_->GetCurrentTrack();
351 mioni_manager_->CreateMultipleIonisedWaterMolecule(
352 eTripleIonisedMolecule, ioni_shell, the_track);
353 // ***************************************************************************
354
355 }
356