| 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 // G4DNATripleIonisationModel.cc
27 //
28 // Created at 2024/04/03 (Thu.)
29 // Author: Shogo OKADA @KEK-CRC (shogo.okada@
30 //
31 // Reference: J.Meesungnoen et. al, DOI: 10.1
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::G4DNATripleIonisat
59 const G4ParticleDefinition* p, const G4Strin
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
69 {
70 if (verbose_level_ > 3) {
71 G4cout << "Calling G4DNATripleIonisationMo
72 }
73
74 proton_def_ = G4Proton::ProtonDefinition();
75 alpha_def_ = G4DNAGenericIonsManager::Insta
76 carbon_def_ = G4IonTable::GetIonTable()->Get
77
78 constexpr G4double kScaleFactor = 1.0 * m *
79
80 mioni_manager_ = new G4DNAMultipleIonisation
81
82 G4double Z{0.0}, A{0.0};
83 G4String alpha_param_file{"dna/multipleionis
84
85 if (particle == proton_def_) {
86
87 // ***************************************
88 // for protons
89 const auto& proton = proton_def_->GetParti
90 elow_tab_[proton] = model_elow_tab_[1];
91 eupp_tab_[proton] = 3.0 * MeV;
92
93 // load cross-section data for single ioni
94 auto xs_proton = new G4DNACrossSectionData
95 new G4LogLogInterpol
96 xs_proton->LoadData("dna/sigma_ionisation_
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/multipleionisati
105 }
106
107 Z = static_cast<G4double>(proton_def_->Get
108 A = static_cast<G4double>(proton_def_->Get
109
110 } else if (particle == alpha_def_) {
111
112 //****************************************
113 // for alpha particles
114 const auto& alpha = alpha_def_->GetParticl
115 elow_tab_[alpha] = model_elow_tab_[4];
116 eupp_tab_[alpha] = 23.0 * MeV;
117
118 // load cross-section data for single ioni
119 auto xs_alpha = new G4DNACrossSectionDataS
120 new G4LogLogInterpolat
121 xs_alpha->LoadData("dna/sigma_ionisation_a
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/multipleionisati
130 }
131
132 Z = static_cast<G4double>(alpha_def_->GetA
133 A = static_cast<G4double>(alpha_def_->GetA
134
135 } else if (particle == G4GenericIon::Generic
136
137 // ***************************************
138 // for carbon ions
139 const auto& carbon = carbon_def_->GetParti
140 elow_tab_[carbon] = model_elow_tab_[5] * c
141 eupp_tab_[carbon] = 120.0 * MeV;
142
143 // load cross-section data for single ioni
144 auto xs_carbon = new G4DNACrossSectionData
145 new G4LogLogInterpolat
146 xs_carbon->LoadData("dna/sigma_ionisation_
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/multipleionisati
155 }
156
157 Z = static_cast<G4double>(carbon_def_->Get
158 A = static_cast<G4double>(carbon_def_->Get
159
160 }
161
162 // load alpha parameter
163 mioni_manager_->LoadAlphaParam(alpha_param_f
164
165 if (verbose_level_ > 0) {
166 G4cout << "G4DNATripleIonisationModel is i
167 << "Energy range: "
168 << LowEnergyLimit() / eV << " eV -
169 << HighEnergyLimit() / keV << " keV
170 << particle->GetParticleName()
171 << G4endl;
172 }
173
174 water_density_ = G4DNAMolecularMaterial::Ins
175 G4Material::GetMaterial(
176
177 atom_deex_ = G4LossTableManager::Instance()-
178
179 if (is_initialized_) { return; }
180
181 particle_change_ = GetParticleChangeForGamma
182 is_initialized_ = true;
183 }
184
185 //--------------------------------------------
186 G4double G4DNATripleIonisationModel::CrossSect
187 const G4Material* material, const G4Particle
188 G4double ekin, G4double, G4double)
189 {
190
191 if (verbose_level_ > 3) {
192 G4cout << "Calling G4DNATripleIonisationMo
193 << G4endl;
194 }
195
196 // Calculate total cross section for model
197
198 if (pdef != proton_def_ && pdef != alpha_def
199 return 0.0;
200 }
201
202 static G4double water_dens = (*water_density
203
204 const auto& pname = pdef->GetParticleName();
205
206 const auto low_energy_lim = GetLowEnergyLimi
207 const auto upp_energy_lim = GetUppEnergyLimi
208
209 G4double sigma{0.0};
210 if (ekin <= upp_energy_lim) {
211
212 if (ekin < low_energy_lim) { ekin = low_en
213
214 CrossSectionDataTable::iterator pos = xs_t
215 if (pos == xs_tab_.end()) {
216 G4Exception("G4DNATripleIonisationModel:
217 "em0002", FatalException,
218 "Model not applicable to par
219 }
220
221 G4DNACrossSectionDataSet* table = pos->sec
222 if (table != nullptr) {
223 const auto a = mioni_manager_->GetAlphaP
224 sigma = table->FindValue(ekin) * a * a;
225 }
226
227 }
228
229 if (verbose_level_ > 2) {
230
231 std::stringstream msg;
232
233 msg << "----------------------------------
234 msg << " G4DNATripleIonisationModel - XS I
235 msg << " - Kinetic energy(eV): " << ekin/
236 << pdef->GetParticleName() << "\n";
237 msg << " - Cross section per water molecu
238 << sigma / cm / cm << "\n";
239 msg << " - Cross section per water molecu
240 << sigma * water_dens / (1.0 / cm) <<
241 msg << " G4DNATripleIonisationModel - XS I
242 msg << "----------------------------------
243
244 G4cout << msg.str() << G4endl;
245
246 }
247
248 return (sigma * water_dens);
249
250 }
251
252 //--------------------------------------------
253 void G4DNATripleIonisationModel::SampleSeconda
254 std::vector<G4DynamicParticle*>* vsec, const
255 const G4DynamicParticle* particle, G4double,
256 {
257
258 if (verbose_level_ > 3) {
259 G4cout << "Calling SampleSecondaries() of
260 << G4endl;
261 }
262
263 // get the definition for this parent partic
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 = GetLowEnergyLimi
274
275 // *****************************************
276 // stop the transportation process of this p
277 // if its kinetic energy is below the lower
278 if (ekin < low_energy_lim) {
279 particle_change_->SetProposedKineticEnergy
280 particle_change_->ProposeTrackStatus(fStop
281 particle_change_->ProposeLocalEnergyDeposi
282 return;
283 }
284 // *****************************************
285
286 constexpr G4int kNumSecondaries = 3;
287 constexpr G4double kDeltaTheta = pi * 0.666
288
289 G4int ioni_shell[kNumSecondaries] = {0, 0, 0
290 G4double shell_energy[kNumSecondaries];
291
292 auto scale_param = mioni_manager_->GetAlphaP
293 scale_param *= scale_param;
294
295 G4bool is_continue{true};
296 while (1) {
297 ioni_shell[0] = RandomSelect(ekin, scale_p
298 ioni_shell[1] = RandomSelect(ekin, scale_p
299 ioni_shell[2] = RandomSelect(ekin, scale_p
300 is_continue = (ioni_shell[0] == ioni_shell
301 ioni_shell[1] == ioni_shell
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.Ionisa
308 tot_ioni_energy += shell_energy[i];
309 }
310
311 if (ekin < tot_ioni_energy || tot_ioni_energ
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, cou
319 theta, ph
320 theta += kDeltaTheta;
321 }
322
323 // This should never happen
324 if (mioni_manager_->CheckShellEnergy(eTriple
325 G4Exception("G4DNATripleIonisatioModel::Sa
326 "em2050", FatalException, "Negative lo
327 }
328
329 // *****************************************
330 // update kinematics for this parent particl
331 const auto primary_dir = particle->GetMoment
332 particle_change_->ProposeMomentumDirection(p
333
334 const auto scattered_energy = ekin - tot_ion
335
336 // update total amount of shell energy
337 tot_ioni_energy = shell_energy[0] + shell_en
338
339 if (stat_code_) {
340 particle_change_->SetProposedKineticEnergy
341 particle_change_->ProposeLocalEnergyDeposi
342 ekin - scatter
343 } else {
344 particle_change_->SetProposedKineticEnergy
345 particle_change_->ProposeLocalEnergyDeposi
346 }
347
348 // *****************************************
349 // generate triple-ionized water molecules (
350 const auto the_track = particle_change_->Get
351 mioni_manager_->CreateMultipleIonisedWaterMo
352 eTripleIonisedMolecule, io
353 // *****************************************
354
355 }
356