| 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 // G4DNAQuadrupleIonisationModel.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 "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::G4DNAQuadrupleI
61 const G4ParticleDefinition* p, const G4Strin
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
71 {
72 if (verbose_level_ > 3) {
73 G4cout << "Calling G4DNAQuadrupleIonisatio
74 }
75
76 proton_def_ = G4Proton::ProtonDefinition();
77 alpha_def_ = G4DNAGenericIonsManager::Insta
78 carbon_def_ = G4IonTable::GetIonTable()->Get
79
80 constexpr G4double kScaleFactor = 1.0 * m *
81
82 mioni_manager_ = new G4DNAMultipleIonisation
83
84 G4double Z{0.0}, A{0.0};
85 G4String alpha_param_file{"dna/multipleionis
86
87 if (particle == proton_def_) {
88
89 // ***************************************
90 // for protons
91 const auto& proton = proton_def_->GetParti
92 elow_tab_[proton] = model_elow_tab_[1];
93 eupp_tab_[proton] = 3.0 * MeV;
94
95 // load cross-section data for single ioni
96 auto xs_proton = new G4DNACrossSectionData
97 new G4LogLogInterpol
98 xs_proton->LoadData("dna/sigma_ionisation_
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/multipleionisati
107 }
108
109 Z = static_cast<G4double>(proton_def_->Get
110 A = static_cast<G4double>(proton_def_->Get
111
112 } else if (particle == alpha_def_) {
113
114 //****************************************
115 // for alpha particles
116 const auto& alpha = alpha_def_->GetParticl
117 elow_tab_[alpha] = model_elow_tab_[4];
118 eupp_tab_[alpha] = 23.0 * MeV;
119
120 // load cross-section data for single ioni
121 auto xs_alpha = new G4DNACrossSectionDataS
122 new G4LogLogInterpolat
123 xs_alpha->LoadData("dna/sigma_ionisation_a
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/multipleionisati
132 }
133
134 Z = static_cast<G4double>(alpha_def_->GetA
135 A = static_cast<G4double>(alpha_def_->GetA
136
137 } else if (particle == G4GenericIon::Generic
138
139 // ***************************************
140 // for carbon ions
141 const auto& carbon = carbon_def_->GetParti
142 elow_tab_[carbon] = model_elow_tab_[5] * c
143 eupp_tab_[carbon] = 120.0 * MeV;
144
145 // load cross-section data for single ioni
146 auto xs_carbon = new G4DNACrossSectionData
147 new G4LogLogInterpolat
148 xs_carbon->LoadData("dna/sigma_ionisation_
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/multipleionisati
157 }
158
159 Z = static_cast<G4double>(carbon_def_->Get
160 A = static_cast<G4double>(carbon_def_->Get
161
162 }
163
164 // load alpha parameter
165 mioni_manager_->LoadAlphaParam(alpha_param_f
166
167 if (verbose_level_ > 0) {
168 G4cout << "G4DNAQuadrupleIonisationModel i
169 << "Energy range: "
170 << LowEnergyLimit() / eV << " eV -
171 << HighEnergyLimit() / keV << " keV
172 << particle->GetParticleName()
173 << G4endl;
174 }
175
176 water_density_ = G4DNAMolecularMaterial::Ins
177 G4Material::GetMaterial(
178
179 atom_deex_ = G4LossTableManager::Instance()-
180
181 if (is_initialized_) { return; }
182
183 particle_change_ = GetParticleChangeForGamma
184 is_initialized_ = true;
185 }
186
187 //--------------------------------------------
188 G4double G4DNAQuadrupleIonisationModel::CrossS
189 const G4Material* material, const G4Particle
190 G4double ekin, G4double, G4double)
191 {
192
193 if (verbose_level_ > 3) {
194 G4cout << "Calling G4DNAQuadrupleIonisatio
195 << G4endl;
196 }
197
198 // Calculate total cross section for model
199
200 if (pdef != proton_def_ && pdef != alpha_def
201 return 0.0;
202 }
203
204 static G4double water_dens = (*water_density
205
206 const auto& pname = pdef->GetParticleName();
207
208 const auto low_energy_lim = GetLowEnergyLimi
209 const auto upp_energy_lim = GetUppEnergyLimi
210
211 G4double sigma{0.0};
212 if (ekin <= upp_energy_lim) {
213
214 if (ekin < low_energy_lim) { ekin = low_en
215
216 CrossSectionDataTable::iterator pos = xs_t
217 if (pos == xs_tab_.end()) {
218 G4Exception("G4DNAQuadrupleIonisationMod
219 "em0002", FatalException,
220 "Model not applicable to par
221 }
222
223 G4DNACrossSectionDataSet* table = pos->sec
224 if (table != nullptr) {
225 auto scale_param = mioni_manager_->GetAl
226 scale_param = ::g4pow->powA(scale_param,
227 sigma = table->FindValue(ekin) * scale_p
228 }
229
230 }
231
232 if (verbose_level_ > 2) {
233
234 std::stringstream msg;
235
236 msg << "----------------------------------
237 msg << " G4DNAQuadrupleIonisationModel - X
238 msg << " - Kinetic energy(eV): " << ekin/
239 << pdef->GetParticleName() << "\n";
240 msg << " - Cross section per water molecu
241 << sigma / cm / cm << "\n";
242 msg << " - Cross section per water molecu
243 << sigma * water_dens / (1.0 / cm) <<
244 msg << " G4DNAQuadrupleIonisationModel - X
245 msg << "----------------------------------
246
247 G4cout << msg.str() << G4endl;
248
249 }
250
251 return (sigma * water_dens);
252
253 }
254
255 //--------------------------------------------
256 void G4DNAQuadrupleIonisationModel::SampleSeco
257 std::vector<G4DynamicParticle*>* vsec, const
258 const G4DynamicParticle* particle, G4double,
259 {
260
261 if (verbose_level_ > 3) {
262 G4cout << "Calling SampleSecondaries() of
263 << G4endl;
264 }
265
266 // get the definition for this parent partic
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 = GetLowEnergyLimi
277
278 // *****************************************
279 // stop the transportation process of this p
280 // if its kinetic energy is below the lower
281 if (ekin < low_energy_lim) {
282 particle_change_->SetProposedKineticEnergy
283 particle_change_->ProposeTrackStatus(fStop
284 particle_change_->ProposeLocalEnergyDeposi
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
293 G4double shell_energy[kNumSecondaries];
294
295 auto scale_param = mioni_manager_->GetAlphaP
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_p
301 ioni_shell[1] = RandomSelect(ekin, scale_p
302 ioni_shell[2] = RandomSelect(ekin, scale_p
303 ioni_shell[3] = RandomSelect(ekin, scale_p
304 is_continue =
305 (ioni_shell[0] == ioni_shell[1] && ioni_
306 (ioni_shell[1] == ioni_shell[2] && ioni_
307 (ioni_shell[2] == ioni_shell[3] && ioni_
308 (ioni_shell[3] == ioni_shell[0] && ioni_
309 (ioni_shell[0] == ioni_shell[1] && ioni_
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.Ionisa
317 tot_ioni_energy += shell_energy[i];
318 }
319
320 if (ekin < tot_ioni_energy || tot_ioni_energ
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, cou
328 theta, ph
329 theta += kDeltaTheta;
330 }
331
332 // This should never happen
333 if (mioni_manager_->CheckShellEnergy(eQuadru
334 shell_e
335 G4Exception("G4DNAQuadrupleIonisatioModel:
336 "em2050", FatalException, "Negative lo
337 }
338
339 // *****************************************
340 // update kinematics for this parent particl
341 const auto primary_dir = particle->GetMoment
342 particle_change_->ProposeMomentumDirection(p
343
344 const auto scattered_energy = ekin - tot_ion
345
346 // update total amount of shell energy
347 tot_ioni_energy = shell_energy[0] + shell_en
348 shell_energy[2] + shell_en
349
350 if (stat_code_) {
351 particle_change_->SetProposedKineticEnergy
352 particle_change_->ProposeLocalEnergyDeposi
353 ekin - scatter
354 } else {
355 particle_change_->SetProposedKineticEnergy
356 particle_change_->ProposeLocalEnergyDeposi
357 }
358
359 // *****************************************
360 // generate triple-ionized water molecules (
361 const auto the_track = particle_change_->Get
362 mioni_manager_->CreateMultipleIonisedWaterMo
363 eQuadrupleIonisedMolecule,
364 // *****************************************
365
366 }
367