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