| 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 // G4DNAMultipleIonisationManager.cc
27 //
28 // Created at 2024/04/03 (Thu.)
29 // Author: Shogo OKADA @KEK-CRC (shogo.okada@
30 //
31
32 #include "G4DNAMultipleIonisationManager.hh"
33 #include "G4FindDataDir.hh"
34 #include "G4Scheduler.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4Molecule.hh"
37 #include "G4VITTrackHolder.hh"
38 #include "G4H2O.hh"
39 #include "G4DNAChemistryManager.hh"
40
41 //--------------------------------------------
42 void G4DNAMultipleIonisationManager::CreateMul
43 MultipleIonisedModification mod, G4int* shel
44 const G4Track* incoming_track)
45 {
46 if (!G4DNAChemistryManager::IsActivated()) {
47
48 G4int num_shells{0};
49 switch (mod) {
50 case eDoubleIonisedMolecule:
51 num_shells = 2;
52 break;
53 case eTripleIonisedMolecule:
54 num_shells = 3;
55 break;
56 case eQuadrupleIonisedMolecule:
57 num_shells = 4;
58 break;
59 default: // never happen
60 return;
61 }
62
63 auto* H2O = new G4Molecule(G4H2O::Definition
64 for (G4int i = 0; i < num_shells; i++) {
65 H2O->IonizeMolecule(4 - shell_level[i]);
66 }
67
68 constexpr G4double kT0 = 1.0 * picosecond;
69 auto* H2O_track = H2O->BuildTrack(kT0, incom
70 H2O_track->SetParentID(incoming_track->GetTr
71 H2O_track->SetTrackStatus(fStopButAlive);
72 H2O_track->SetKineticEnergy(0.0);
73
74 G4VITTrackHolder::Instance()->Push(H2O_track
75 }
76
77 //--------------------------------------------
78 G4bool G4DNAMultipleIonisationManager::CheckSh
79 MultipleIonisedModification mod, G4double*
80 {
81 G4int num_shells{0};
82 switch (mod) {
83 case eDoubleIonisedMolecule:
84 num_shells = 2;
85 break;
86 case eTripleIonisedMolecule:
87 num_shells = 3;
88 break;
89 case eQuadrupleIonisedMolecule:
90 num_shells = 4;
91 break;
92 default: // never happen
93 break;
94 }
95
96 G4bool stop_process{false};
97 for (G4int i = 0; i < num_shells; i++) {
98 if (shell_energy[i] < 0.0) {
99 stop_process = true;
100 break;
101 }
102 }
103
104 return stop_process;
105 }
106
107 //--------------------------------------------
108 void G4DNAMultipleIonisationManager::LoadAlpha
109 const G4String& filepath, G4double Z, G4doub
110 {
111 const char* path = G4FindDataDir("G4LEDATA")
112 if (path == nullptr) {
113 G4Exception("G4DNAMultipleIonisationManage
114 FatalException,"G4LEDATA envir
115 }
116
117 std::stringstream fullpath;
118 fullpath << path << "/" << filepath;
119
120 std::fstream fin(fullpath.str());
121
122 G4double e, a;
123 std::string line = "";
124 while (getline(fin, line)) {
125 std::stringstream ss;
126 ss << line;
127 ss >> e >> a;
128 Eion_.push_back(e * Z * A * MeV);
129 alpha_.push_back(a);
130 }
131
132 num_node_ = (G4int)Eion_.size();
133 fin.close();
134 }
135
136 //--------------------------------------------
137 G4double G4DNAMultipleIonisationManager::GetAl
138 {
139 auto find_lower_bound = [this](G4double e) {
140 auto low = 0;
141 auto upp = num_node_ - 1;
142 if (e < Eion_[0]) { return low; }
143 while (low <= upp) {
144 const auto mid = static_cast<int>((low +
145 if (e < Eion_[mid]) { upp = mid - 1; }
146 else { low = mid + 1; }
147 if (upp < 0) { upp = 0; }
148 }
149 return upp;
150 };
151
152 auto interp_log_log = [this](G4int bin1, G4d
153 if (e < Eion_[0]) { return alpha_[0]; }
154 const auto num_bin = num_node_ - 1;
155 const auto bin2 = bin1 + 1;
156 G4double value{0.0};
157 if (bin2 <= num_bin) {
158 auto log10_e = std::log10(e);
159 auto log10_e1 = Eion_[bin1];
160 auto log10_e2 = Eion_[bin2];
161 auto log10_a1 = alpha_[bin1];
162 auto log10_a2 = alpha_[bin2];
163 if (log10_a1 != 0.0 && log10_a2 != 0.0)
164 log10_e1 = std::log10(log10_e1);
165 log10_e2 = std::log10(log10_e2);
166 log10_a1 = std::log10(log10_a1);
167 log10_a2 = std::log10(log10_a2);
168 value = log10_a1 + (log10_a2 - log10_a
169 * (log10_e - log10_e1) / (lo
170 value = std::pow(10.0, value);
171 }
172 } else {
173 value = alpha_[num_bin];
174 }
175 return value;
176 };
177
178 const auto bin1 = find_lower_bound(energy);
179 return interp_log_log(bin1, energy);
180 }
181