| 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 // CPA100 elastic model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C.
29 //
30 // Users are requested to cite the following p
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Do
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terr
33 // M. Bardies, N. Lampe, S. Incerti, Phys. M
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bor
37 //
38 // 15.01.2014: creation
39 //
40 // Based on the study by S. Zein et. al. Nucl.
41 // 1/2/2023 : Hoang added modification
42
43 #include "G4DNACPA100ElasticModel.hh"
44
45 #include "G4DNAMaterialManager.hh"
46 #include "G4DNAMolecularMaterial.hh"
47 #include "G4EnvironmentUtils.hh"
48 #include "G4SystemOfUnits.hh"
49 using namespace std;
50 G4DNACPA100ElasticModel::G4DNACPA100ElasticMod
51 : G4VDNAModel(nam, "all")
52 {
53 fpGuanine = G4Material::GetMaterial("G4_GUAN
54 fpG4_WATER = G4Material::GetMaterial("G4_WAT
55 fpDeoxyribose = G4Material::GetMaterial("G4_
56 fpCytosine = G4Material::GetMaterial("G4_CYT
57 fpThymine = G4Material::GetMaterial("G4_THYM
58 fpAdenine = G4Material::GetMaterial("G4_ADEN
59 fpPhosphate = G4Material::GetMaterial("G4_PH
60 fpParticle = G4Electron::ElectronDefinition(
61 }
62
63 void G4DNACPA100ElasticModel::Initialise(const
64 const
65 {
66 if (isInitialised) {
67 return;
68 }
69 if (verboseLevel > 3) {
70 G4cout << "Calling G4DNACPA100ExcitationMo
71 }
72
73 if (!G4DNAMaterialManager::Instance()->IsLoc
74 if (p != fpParticle) {
75 std::ostringstream oss;
76 oss << " Model is not applied for this p
77 G4Exception("G4DNACPA100ElasticModel::G4
78 oss.str().c_str());
79 }
80
81 const char* path = G4FindDataDir("G4LEDATA
82
83 if (path == nullptr) {
84 G4Exception("G4DNACPA100ElasticModel::In
85 "G4LEDATA environment variab
86 return;
87 }
88
89 std::size_t index;
90 if (fpG4_WATER != nullptr) {
91 index = fpG4_WATER->GetIndex();
92 fLevels[index] = 1.214e-4;
93 AddCrossSectionData(index, p, "dna/sigma
94 "dna/sigmadiff_cumul
95 SetLowELimit(index, p, 11. * eV);
96 SetHighELimit(index, p, 255955. * eV);
97 }
98 if (fpGuanine != nullptr) {
99 index = fpGuanine->GetIndex();
100 fLevels[index] = 1.4504480e-05;
101 AddCrossSectionData(index, p, "dna/sigma
102 "dna/sigmadiff_cumul
103 SetLowELimit(index, p, 11 * eV);
104 SetHighELimit(index, p, 1 * MeV);
105 }
106 if (fpDeoxyribose != nullptr) {
107 index = fpDeoxyribose->GetIndex();
108 fLevels[index] = 1.6343100e-05;
109 AddCrossSectionData(index, p, "dna/sigma
110 "dna/sigmadiff_cumul
111 SetLowELimit(index, p, 11 * eV);
112 SetHighELimit(index, p, 1 * MeV);
113 }
114 if (fpCytosine != nullptr) {
115 index = fpCytosine->GetIndex();
116 fLevels[index] = 1.9729660e-05;
117 AddCrossSectionData(index, p, "dna/sigma
118 "dna/sigmadiff_cumul
119 SetLowELimit(index, p, 11 * eV);
120 SetHighELimit(index, p, 1 * MeV);
121 }
122 if (fpThymine != nullptr) {
123 index = fpThymine->GetIndex();
124 fLevels[index] = 1.7381300e-05;
125 AddCrossSectionData(index, p, "dna/sigma
126 "dna/sigmadiff_cumul
127 SetLowELimit(index, p, 11 * eV);
128 SetHighELimit(index, p, 1 * MeV);
129 }
130 if (fpAdenine != nullptr) {
131 index = fpAdenine->GetIndex();
132 fLevels[index] = 1.6221800e-05;
133 AddCrossSectionData(index, p, "dna/sigma
134 "dna/sigmadiff_cumul
135 SetLowELimit(index, p, 11 * eV);
136 SetHighELimit(index, p, 1 * MeV);
137 }
138 if (fpPhosphate != nullptr) {
139 index = fpPhosphate->GetIndex();
140 fLevels[index] = 2.2369600e-05;
141 AddCrossSectionData(index, p, "dna/sigma
142 "dna/sigmadiff_cumul
143 SetLowELimit(index, p, 11 * eV);
144 SetHighELimit(index, p, 1 * MeV);
145 }
146
147 // Load data
148 LoadCrossSectionData(p);
149 G4DNAMaterialManager::Instance()->SetMaste
150 fpModelData = this;
151 }
152 else {
153 auto dataModel = dynamic_cast<G4DNACPA100E
154 G4DNAMaterialManager::Instance()->GetMod
155 if (dataModel == nullptr) {
156 G4cout << "G4DNACPA100ElasticModel::Cros
157 G4Exception("G4DNACPA100ElasticModel::Cr
158 "no modelData is registered"
159 }
160 else {
161 fpModelData = dataModel;
162 }
163 }
164
165 fParticleChangeForGamma = GetParticleChangeF
166 isInitialised = true;
167 }
168
169 //....oooOO0OOooo........oooOO0OOooo........oo
170
171 G4double G4DNACPA100ElasticModel::CrossSection
172
173
174 {
175 // Get the name of the current particle
176 const G4String& particleName = p->GetParticl
177 auto materialID = pMaterial->GetIndex();
178
179 // set killBelowEnergy value for current mat
180 fKillBelowEnergy = fpModelData->GetLowELimit
181
182 G4double sigma = 0.;
183
184 if (ekin < fpModelData->GetHighELimit(materi
185 if (ekin < fKillBelowEnergy) {
186 return DBL_MAX;
187 }
188
189 auto tableData = fpModelData->GetData();
190
191 if ((*tableData)[materialID][p] == nullptr
192 G4Exception("G4DNACPA100ElasticModel::Cr
193 "No model is registered");
194 }
195 sigma = (*tableData)[materialID][p]->FindV
196 }
197
198 if (verboseLevel > 2) {
199 auto MolDensity =
200 (*G4DNAMolecularMaterial::Instance()->Ge
201 G4cout << "_______________________________
202 G4cout << "°°° G4DNACPA100ElasticModel
203 G4cout << "°°° Kinetic energy(eV)=" <<
204 G4cout << "°°° lowLim (eV) = " << GetLo
205 << " highLim (eV) : " << GetHighELi
206 G4cout << "°°° Materials = " << (*G4Mat
207 << G4endl;
208 G4cout << "°°° Cross section per molecu
209 G4cout << "°°° Cross section per Phosph
210 << G4endl;
211 G4cout << "°°° G4DNACPA100ElasticModel
212 }
213
214 auto MolDensity =
215 (*G4DNAMolecularMaterial::Instance()->GetN
216 return sigma * MolDensity;
217 }
218
219 //....oooOO0OOooo........oooOO0OOooo........oo
220
221 void G4DNACPA100ElasticModel::SampleSecondarie
222
223
224
225 {
226 G4double electronEnergy0 = aDynamicElectron-
227 auto materialID = couple->GetMaterial()->Get
228 auto p = aDynamicElectron->GetParticleDefini
229
230 if (p != fpParticle) {
231 G4Exception("G4DNACPA100ElasticModel::Samp
232 "This particle is not applied
233 }
234 if (electronEnergy0 < fKillBelowEnergy) {
235 return;
236 }
237 G4double cosTheta = fpModelData->RandomizeCo
238 G4double phi = 2. * CLHEP::pi * G4UniformRan
239
240 const G4ThreeVector& zVers = aDynamicElectro
241
242 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2,
243 G4double sinTheta = std::sqrt(1 - cosTheta *
244
245 CT1 = zVers.z();
246 ST1 = std::sqrt(1. - CT1 * CT1);
247
248 if (ST1 != 0)
249 CF1 = zVers.x() / ST1;
250 else
251 CF1 = std::cos(2. * CLHEP::pi * G4UniformR
252 if (ST1 != 0)
253 SF1 = zVers.y() / ST1;
254 else
255 SF1 = std::sqrt(1. - CF1 * CF1);
256
257 G4double A3, A4, A5, A2, A1;
258
259 A3 = sinTheta * std::cos(phi);
260 A4 = A3 * CT1 + ST1 * cosTheta;
261 A5 = sinTheta * std::sin(phi);
262 A2 = A4 * SF1 + A5 * CF1;
263 A1 = A4 * CF1 - A5 * SF1;
264
265 CT2 = CT1 * cosTheta - ST1 * A3;
266 ST2 = std::sqrt(1. - CT2 * CT2);
267
268 if (ST2 == 0) ST2 = 1E-6;
269 CF2 = A1 / ST2;
270 SF2 = A2 / ST2;
271 G4ThreeVector zPrimeVers(ST2 * CF2, ST2 * SF
272
273 fParticleChangeForGamma->ProposeMomentumDire
274
275 auto EnergyDeposit = fpModelData->GetElastic
276 fParticleChangeForGamma->ProposeLocalEnergyD
277 if (statCode) {
278 fParticleChangeForGamma->SetProposedKineti
279 }
280 else {
281 auto newEnergy = electronEnergy0 - EnergyD
282 fParticleChangeForGamma->SetProposedKineti
283 }
284 }
285
286 G4double G4DNACPA100ElasticModel::Theta(const
287 G4doub
288 {
289 G4double theta, valueT1, valueT2, valueE21,
290 G4double xs11 = 0;
291 G4double xs12 = 0;
292 G4double xs21 = 0;
293 G4double xs22 = 0;
294 if (p == G4Electron::ElectronDefinition()) {
295 if (k == tValuesVec[materialID][p].back())
296 k = k * (1. - 1e-12);
297 }
298 auto t2 =
299 std::upper_bound(tValuesVec[materialID][
300 auto t1 = t2 - 1;
301
302 auto e12 = std::upper_bound(eValuesVect[ma
303 eValuesVect[ma
304 auto e11 = e12 - 1;
305
306 auto e22 = std::upper_bound(eValuesVect[ma
307 eValuesVect[ma
308 auto e21 = e22 - 1;
309
310 valueT1 = *t1;
311 valueT2 = *t2;
312 valueE21 = *e21;
313 valueE22 = *e22;
314 valueE12 = *e12;
315 valueE11 = *e11;
316
317 xs11 = diffCrossSectionData[materialID][p]
318 xs12 = diffCrossSectionData[materialID][p]
319 xs21 = diffCrossSectionData[materialID][p]
320 xs22 = diffCrossSectionData[materialID][p]
321 }
322
323 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x
324 return (0.);
325 }
326
327 theta = QuadInterpolator(valueE11, valueE12,
328 valueT2, k, integrD
329
330 return theta;
331 }
332
333 //....oooOO0OOooo........oooOO0OOooo........oo
334
335 G4double G4DNACPA100ElasticModel::LinLogInterp
336
337 {
338 G4double d1 = std::log(xs1);
339 G4double d2 = std::log(xs2);
340 G4double value = std::exp(d1 + (d2 - d1) * (
341 return value;
342 }
343
344 //....oooOO0OOooo........oooOO0OOooo........oo
345
346 G4double G4DNACPA100ElasticModel::LinLinInterp
347
348 {
349 G4double d1 = xs1;
350 G4double d2 = xs2;
351 G4double value = (d1 + (d2 - d1) * (e - e1)
352 return value;
353 }
354
355 //....oooOO0OOooo........oooOO0OOooo........oo
356
357 G4double G4DNACPA100ElasticModel::LogLogInterp
358
359 {
360 G4double a = (std::log10(xs2) - std::log10(x
361 G4double b = std::log10(xs2) - a * std::log1
362 G4double sigma = a * std::log10(e) + b;
363 G4double value = (std::pow(10., sigma));
364 return value;
365 }
366
367 //....oooOO0OOooo........oooOO0OOooo........oo
368
369 G4double G4DNACPA100ElasticModel::QuadInterpol
370
371
372
373 {
374 // Log-Log
375 /*
376 G4double interpolatedvalue1 = LogLogInterp
377 G4double interpolatedvalue2 = LogLogInterp
378 G4double value = LogLogInterpolate(t1, t2,
379
380
381 // Lin-Log
382 G4double interpolatedvalue1 = LinLogInterp
383 G4double interpolatedvalue2 = LinLogInterp
384 G4double value = LinLogInterpolate(t1, t2,
385 */
386
387 // Lin-Lin
388 G4double interpolatedvalue1 = LinLinInterpol
389 G4double interpolatedvalue2 = LinLinInterpol
390 G4double value = LinLinInterpolate(t1, t2, t
391
392 return value;
393 }
394
395 //....oooOO0OOooo........oooOO0OOooo........oo
396
397 G4double G4DNACPA100ElasticModel::RandomizeCos
398 {
399 G4double integrdiff = 0; // PROBABILITY bet
400 G4double uniformRand = G4UniformRand();
401 integrdiff = uniformRand;
402 G4double cosTheta = 0.;
403 cosTheta = 1 - Theta(G4Electron::ElectronDef
404 // cosTheta = std::cos(theta * CLHEP::pi / 1
405 return cosTheta;
406 }
407 //....oooOO0OOooo........oooOO0OOooo........oo
408
409 void G4DNACPA100ElasticModel::ReadDiffCSFile(c
410 c
411 c
412 {
413 const char* path = G4FindDataDir("G4LEDATA")
414 if (path == nullptr) {
415 G4Exception("G4DNACPA100ElasticModel::Read
416 "G4LEDATA environment variable
417 return;
418 }
419
420 std::ostringstream fullFileName;
421 fullFileName << path << "/" << file << ".dat
422
423 std::ifstream diffCrossSection(fullFileName.
424 // error if file is not there
425 std::stringstream endPath;
426 if (!diffCrossSection) {
427 endPath << "Missing data file: " << file;
428 G4Exception("G4DNACPA100ElasticModel::Init
429 endPath.str().c_str());
430 }
431
432 tValuesVec[materialName][particleName].push_
433
434 G4String line;
435 while (std::getline(diffCrossSection, line))
436 //
437 std::istringstream testIss(line);
438 G4String test;
439 testIss >> test;
440 if (test == "#") {
441 continue;
442 }
443 // check if line is empty
444 if (line.empty()) {
445 continue;
446 }
447 std::istringstream iss(line);
448
449 G4double tDummy;
450 G4double eDummy;
451
452 iss >> tDummy >> eDummy;
453
454 if (tDummy != tValuesVec[materialName][par
455 // Add the current T value
456 tValuesVec[materialName][particleName].p
457 // Make it correspond to a default zero
458 eValuesVect[materialName][particleName][
459 }
460 iss >> diffCrossSectionData[materialName][
461
462 if (eDummy != eValuesVect[materialName][pa
463 eValuesVect[materialName][particleName][
464 }
465 }
466 }
467