| 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 // Authors: S. Meylan and C. Villagrasa (IRSN,
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-
29 //
30
31 #include "G4DNAPTBExcitationModel.hh"
32
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMaterialManager.hh"
35 #include "G4DNAMolecularMaterial.hh"
36 #include "G4SystemOfUnits.hh"
37
38 G4DNAPTBExcitationModel::G4DNAPTBExcitationMod
39
40 : G4VDNAModel(nam, applyToMaterial)
41 {
42 fpTHF = G4Material::GetMaterial("THF", false
43 fpPY = G4Material::GetMaterial("PY", false);
44 fpPU = G4Material::GetMaterial("PU", false);
45 fpTMP = G4Material::GetMaterial("TMP", false
46 fpG4_WATER = G4Material::GetMaterial("G4_WAT
47 fpBackbone_THF = G4Material::GetMaterial("ba
48 fpCytosine_PY = G4Material::GetMaterial("cyt
49 fpThymine_PY = G4Material::GetMaterial("thym
50 fpAdenine_PU = G4Material::GetMaterial("aden
51 fpBackbone_TMP = G4Material::GetMaterial("ba
52 fpGuanine_PU = G4Material::GetMaterial("guan
53 fpN2 = G4Material::GetMaterial("N2", false);
54 // initialisation of mean energy loss for ea
55
56 if (fpTHF != nullptr) {
57 fTableMeanEnergyPTB[fpTHF->GetIndex()] = 8
58 }
59
60 if (fpPY != nullptr) {
61 fTableMeanEnergyPTB[fpPY->GetIndex()] = 7.
62 }
63
64 if (fpPU != nullptr) {
65 fTableMeanEnergyPTB[fpPU->GetIndex()] = 7.
66 }
67 if (fpTMP != nullptr) {
68 fTableMeanEnergyPTB[fpTMP->GetIndex()] = 8
69 }
70
71 if (verboseLevel > 0) {
72 G4cout << "PTB excitation model is constru
73 }
74 }
75
76 //....oooOO0OOooo........oooOO0OOooo........oo
77
78 void G4DNAPTBExcitationModel::Initialise(const
79 const
80 {
81 if (isInitialised) {
82 return;
83 }
84 if (verboseLevel > 3)
85 {
86 G4cout << "Calling G4DNAPTBExcitationModel
87 }
88
89 if (particle != G4Electron::ElectronDefiniti
90 std::ostringstream oss;
91 oss << " Model is not applied for this par
92 G4Exception("G4DNAPTBExcitationModel::Init
93 }
94
95 G4double scaleFactor = 1e-16 * cm * cm;
96 G4double scaleFactorBorn = (1.e-22 / 3.343)
97
98 //******************************************
99 // Cross section data
100 //******************************************
101 std::size_t index;
102 if (fpTHF != nullptr) {
103 index = fpTHF->GetIndex();
104 AddCrossSectionData(index, particle, "dna/
105 SetLowELimit(index, particle, 9. * eV);
106 SetHighELimit(index, particle, 1. * keV);
107 }
108 if (fpPY != nullptr) {
109 index = fpPY->GetIndex();
110 AddCrossSectionData(index, particle, "dna/
111 SetLowELimit(index, particle, 9. * eV);
112 SetHighELimit(index, particle, 1. * keV);
113 }
114
115 if (fpPU != nullptr) {
116 index = fpPU->GetIndex();
117 AddCrossSectionData(index, particle, "dna/
118 SetLowELimit(index, particle, 9. * eV);
119 SetHighELimit(index, particle, 1. * keV);
120 }
121
122 if (fpTMP != nullptr) {
123 index = fpTMP->GetIndex();
124 AddCrossSectionData(index, particle, "dna/
125 SetLowELimit(index, particle, 9. * eV);
126 SetHighELimit(index, particle, 1. * keV);
127 }
128 if (fpG4_WATER != nullptr) {
129 index = fpG4_WATER->GetIndex();
130 AddCrossSectionData(index, particle, "dna/
131 SetLowELimit(index, particle, 9. * eV);
132 SetHighELimit(index, particle, 1. * keV);
133 }
134 // DNA materials
135 //
136 if (fpBackbone_THF != nullptr) {
137 index = fpBackbone_THF->GetIndex();
138 AddCrossSectionData(index, particle, "dna/
139 SetLowELimit(index, particle, 9. * eV);
140 SetHighELimit(index, particle, 1. * keV);
141 }
142 if (fpCytosine_PY != nullptr) {
143 index = fpCytosine_PY->GetIndex();
144 AddCrossSectionData(index, particle, "dna/
145 SetLowELimit(index, particle, 9. * eV);
146 SetHighELimit(index, particle, 1. * keV);
147 }
148 if (fpThymine_PY != nullptr) {
149 index = fpThymine_PY->GetIndex();
150 AddCrossSectionData(index, particle, "dna/
151 SetLowELimit(index, particle, 9. * eV);
152 SetHighELimit(index, particle, 1. * keV);
153 }
154 if (fpAdenine_PU != nullptr) {
155 index = fpAdenine_PU->GetIndex();
156 AddCrossSectionData(index, particle, "dna/
157 SetLowELimit(index, particle, 9. * eV);
158 SetHighELimit(index, particle, 1. * keV);
159 }
160 if (fpGuanine_PU != nullptr) {
161 index = fpGuanine_PU->GetIndex();
162 AddCrossSectionData(index, particle, "dna/
163 SetLowELimit(index, particle, 9. * eV);
164 SetHighELimit(index, particle, 1. * keV);
165 }
166 if (fpBackbone_TMP != nullptr) {
167 index = fpBackbone_TMP->GetIndex();
168 AddCrossSectionData(index, particle, "dna/
169 SetLowELimit(index, particle, 9. * eV);
170 SetHighELimit(index, particle, 1. * keV);
171 }
172 // MPietrzak, adding paths for N2
173 if (fpN2 != nullptr) {
174 index = fpN2->GetIndex();
175 AddCrossSectionData(index, particle, "dna/
176 SetLowELimit(index, particle, 13. * eV);
177 SetHighELimit(index, particle, 1.02 * MeV)
178 }
179 if (!G4DNAMaterialManager::Instance()->IsLoc
180 // Load data
181 LoadCrossSectionData(particle);
182 G4DNAMaterialManager::Instance()->SetMaste
183 fpModelData = this;
184 }
185 else {
186 auto dataModel = dynamic_cast<G4DNAPTBExci
187 G4DNAMaterialManager::Instance()->GetMod
188 if (dataModel == nullptr) {
189 G4cout << "G4DNAPTBExcitationModel::Init
190 G4Exception("G4DNAPTBExcitationModel::In
191 "not good modelData");
192 }
193 else {
194 fpModelData = dataModel;
195 }
196 }
197
198 fParticleChangeForGamma = GetParticleChangeF
199 isInitialised = true;
200 }
201
202 //....oooOO0OOooo........oooOO0OOooo........oo
203
204 G4double G4DNAPTBExcitationModel::CrossSection
205
206
207
208 {
209 // Get the name of the current particle
210 G4String particleName = p->GetParticleName()
211 const std::size_t& MatID = material->GetInde
212 // initialise variables
213 G4double lowLim;
214 G4double highLim;
215 G4double sigma = 0;
216
217 // Get the low energy limit for the current
218 lowLim = fpModelData->GetLowELimit(MatID, p)
219
220 // Get the high energy limit for the current
221 highLim = fpModelData->GetHighELimit(MatID,
222
223 // Check that we are in the correct energy r
224 if (ekin >= lowLim && ekin < highLim) {
225 // Get the map with all the data tables
226 auto Data = fpModelData->GetData();
227
228 if ((*Data)[MatID][p] == nullptr) {
229 G4Exception("G4DNAPTBExcitationModel::Cr
230 "No model is registered");
231 }
232 // Retrieve the cross section value
233 sigma = (*Data)[MatID][p]->FindValue(ekin)
234
235 if (verboseLevel > 2) {
236 G4cout << "_____________________________
237 G4cout << "°°° G4DNAPTBExcitationMode
238 G4cout << "°°° Kinetic energy(eV)=" <
239 G4cout << "°°° Cross section per " <<
240 << G4endl;
241 G4cout << "°°° G4DNAPTBExcitationMode
242 }
243 }
244
245 // Return the cross section value
246 auto MolDensity =
247 (*G4DNAMolecularMaterial::Instance()->GetN
248 return sigma * MolDensity;
249 }
250
251 //....oooOO0OOooo........oooOO0OOooo........oo
252
253 void G4DNAPTBExcitationModel::SampleSecondarie
254
255
256
257 {
258 const std::size_t& materialID = (std::size_t
259
260 // Get the incident particle kinetic energy
261 G4double k = aDynamicParticle->GetKineticEne
262 // Get the particle name
263 const auto& particle = aDynamicParticle->Get
264 // Get the energy limits
265 G4double lowLim = fpModelData->GetLowELimit(
266 G4double highLim = fpModelData->GetHighELimi
267
268 // Check if we are in the correct energy ran
269 if (k >= lowLim && k < highLim) {
270 if (fpN2 != nullptr && materialID == fpN2-
271 // Retrieve the excitation energy for th
272 G4int level = fpModelData->RandomSelectS
273 G4double excitationEnergy = ptbExcitatio
274
275 // Calculate the new energy of the parti
276 G4double newEnergy = k - excitationEnerg
277
278 // Check that the new energy is above ze
279 // Otherwise, do nothing.
280 if (newEnergy > 0) {
281 fParticleChangeForGamma->ProposeMoment
282 fParticleChangeForGamma->SetProposedKi
283 fParticleChangeForGamma->ProposeLocalE
284 G4double ioniThres = ptbIonisationStru
285 // if excitation energy greater than i
286 if ((excitationEnergy > ioniThres) &&
287 fParticleChangeForGamma->ProposeLoca
288 // energy of ejected electron
289 G4double secondaryKinetic = excitati
290 // random direction
291 G4double cosTheta = 2 * G4UniformRan
292 G4double sinTheta = std::sqrt(1. - c
293 G4double ux = sinTheta * std::cos(ph
294 G4ThreeVector deltaDirection(ux, uy,
295 // Create the new particle with its
296 auto dp = new G4DynamicParticle(G4El
297 fvect->push_back(dp);
298 }
299 }
300 else {
301 G4ExceptionDescription description;
302 description << "Kinetic energy <= 0 at
303 G4Exception("G4DNAPTBExcitationModel::
304 }
305 }
306 else if (fpG4_WATER == nullptr || material
307 // Retrieve the excitation energy for th
308 G4double excitationEnergy = fTableMeanEn
309 // Calculate the new energy of the parti
310 G4double newEnergy = k - excitationEnerg
311 // Check that the new energy is above ze
312 // Otherwise, do nothing.
313 if (newEnergy > 0) {
314 fParticleChangeForGamma->ProposeMoment
315 fParticleChangeForGamma->SetProposedKi
316 fParticleChangeForGamma->ProposeLocalE
317 }
318 else {
319 G4ExceptionDescription description;
320 description << "Kinetic energy <= 0 at
321 G4Exception("G4DNAPTBExcitationModel::
322 }
323 }
324 else {
325 G4int level = RandomSelectShell(k, parti
326 G4double excitationEnergy = waterStructu
327 G4double newEnergy = k - excitationEnerg
328
329 if (newEnergy > 0) {
330 fParticleChangeForGamma->ProposeMoment
331 fParticleChangeForGamma->SetProposedKi
332 fParticleChangeForGamma->ProposeLocalE
333 const G4Track* theIncomingTrack = fPar
334 G4DNAChemistryManager::Instance()->Cre
335
336 }
337 else {
338 G4ExceptionDescription description;
339 description << "Kinetic energy <= 0 at
340 G4Exception("G4DNAPTBExcitationModel::
341 }
342 }
343 }
344 }
345