| 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 "G4DNAPTBElasticModel.hh"
32
33 #include "G4DNAChampionElasticModel.hh"
34 #include "G4DNAMaterialManager.hh"
35 #include "G4DNAMolecularMaterial.hh"
36 #include "G4Proton.hh"
37 #include "G4SystemOfUnits.hh"
38
39 G4DNAPTBElasticModel::G4DNAPTBElasticModel(con
40 con
41 : G4VDNAModel(nam, applyToMaterial)
42 {
43 if (verboseLevel > 0) {
44 G4cout << "PTB Elastic model is constructe
45 }
46 fpTHF = G4Material::GetMaterial("THF", false
47 fpPY = G4Material::GetMaterial("PY", false);
48 fpPU = G4Material::GetMaterial("PU", false);
49 fpTMP = G4Material::GetMaterial("TMP", false
50 fpG4_WATER = G4Material::GetMaterial("G4_WAT
51 fpBackbone_THF = G4Material::GetMaterial("ba
52 fpCytosine_PY = G4Material::GetMaterial("cyt
53 fpThymine_PY = G4Material::GetMaterial("thym
54 fpAdenine_PU = G4Material::GetMaterial("aden
55 fpBackbone_TMP = G4Material::GetMaterial("ba
56 fpGuanine_PU = G4Material::GetMaterial("guan
57 fpN2 = G4Material::GetMaterial("N2", false);
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oo
61
62 void G4DNAPTBElasticModel::Initialise(const G4
63 const G4
64 {
65 if (isInitialised) {
66 return;
67 }
68 if (verboseLevel > 3)
69 {
70 G4cout << "Calling G4DNAPTBElasticModel::I
71 }
72 if (particle != G4Electron::ElectronDefiniti
73 std::ostringstream oss;
74 oss << " Model is not applied for this par
75 G4Exception("G4DNAPTBElasticModel::G4DNAPT
76 oss.str().c_str());
77 }
78 G4double scaleFactor = 1e-16 * cm * cm;
79 //******************************************
80 // Cross section data
81 //******************************************
82
83 std::size_t index;
84 // MPietrzak, adding paths for N2
85 if (fpN2 != nullptr) {
86 index = fpN2->GetIndex();
87 AddCrossSectionData(index, particle, "dna/
88 "dna/sigmadiff_cumulat
89 SetLowELimit(index, particle, 10 * eV);
90 SetHighELimit(index, particle, 1.02 * MeV)
91 }
92 // MPietrzak
93
94 if (fpTHF != nullptr) {
95 index = fpTHF->GetIndex();
96 AddCrossSectionData(index, particle, "dna/
97 "dna/sigmadiff_cumulat
98 SetLowELimit(index, particle, 10 * eV);
99 SetHighELimit(index, particle, 1 * keV);
100 }
101
102 if (fpPY != nullptr) {
103 index = fpPY->GetIndex();
104 AddCrossSectionData(index, particle, "dna/
105 "dna/sigmadiff_cumulat
106 SetLowELimit(index, particle, 10 * eV);
107 SetHighELimit(index, particle, 1 * keV);
108 }
109
110 if (fpPU != nullptr) {
111 index = fpPU->GetIndex();
112 AddCrossSectionData(index, particle, "dna/
113 "dna/sigmadiff_cumulat
114 SetLowELimit(index, particle, 10 * eV);
115 SetHighELimit(index, particle, 1 * keV);
116 }
117
118 if (fpTMP != nullptr) {
119 index = fpTMP->GetIndex();
120 AddCrossSectionData(index, particle, "dna/
121 "dna/sigmadiff_cumulat
122 SetLowELimit(index, particle, 10 * eV);
123 SetHighELimit(index, particle, 1 * keV);
124 }
125 //????
126 if (fpG4_WATER != nullptr) {
127 index = fpG4_WATER->GetIndex();
128 AddCrossSectionData(index, particle, "dna/
129 "dna/sigmadiff_cumulat
130 SetLowELimit(index, particle, 10 * eV);
131 SetHighELimit(index, particle, 1 * keV);
132 }
133 // DNA materials
134 //
135 if (fpBackbone_THF != nullptr) {
136 index = fpBackbone_THF->GetIndex();
137 AddCrossSectionData(index, particle, "dna/
138 "dna/sigmadiff_cumulat
139 SetLowELimit(index, particle, 10 * eV);
140 SetHighELimit(index, particle, 1 * keV);
141 }
142
143 if (fpCytosine_PY != nullptr) {
144 index = fpCytosine_PY->GetIndex();
145 AddCrossSectionData(index, particle, "dna/
146 "dna/sigmadiff_cumulat
147 SetLowELimit(index, particle, 10 * eV);
148 SetHighELimit(index, particle, 1 * keV);
149 }
150
151 if (fpThymine_PY != nullptr) {
152 index = fpThymine_PY->GetIndex();
153 AddCrossSectionData(index, particle, "dna/
154 "dna/sigmadiff_cumulat
155 SetLowELimit(index, particle, 10 * eV);
156 SetHighELimit(index, particle, 1 * keV);
157 }
158
159 if (fpAdenine_PU != nullptr) {
160 index = fpAdenine_PU->GetIndex();
161 AddCrossSectionData(index, particle, "dna/
162 "dna/sigmadiff_cumulat
163 SetLowELimit(index, particle, 10 * eV);
164 SetHighELimit(index, particle, 1 * keV);
165 }
166 if (fpGuanine_PU != nullptr) {
167 index = fpGuanine_PU->GetIndex();
168 AddCrossSectionData(index, particle, "dna/
169 "dna/sigmadiff_cumulat
170 SetLowELimit(index, particle, 10 * eV);
171 SetHighELimit(index, particle, 1 * keV);
172 }
173
174 if (fpBackbone_TMP != nullptr) {
175 index = fpBackbone_TMP->GetIndex();
176 AddCrossSectionData(index, particle, "dna/
177 "dna/sigmadiff_cumulat
178 SetLowELimit(index, particle, 10 * eV);
179 SetHighELimit(index, particle, 1 * keV);
180 }
181
182 if (!G4DNAMaterialManager::Instance()->IsLoc
183 // Load the data
184 LoadCrossSectionData(particle);
185 G4DNAMaterialManager::Instance()->SetMaste
186 fpModelData = this;
187 }
188 else {
189 auto dataModel = dynamic_cast<G4DNAPTBElas
190 G4DNAMaterialManager::Instance()->GetMod
191 if (dataModel == nullptr) {
192 G4cout << "G4DNAPTBElasticModel::Initial
193 G4Exception("G4DNAPTBElasticModel::Initi
194 "not good modelData");
195 }
196 else {
197 fpModelData = dataModel;
198 }
199 }
200
201 if (verboseLevel > 2) {
202 G4cout << "Loaded cross section files for
203 }
204
205 fParticleChangeForGamma = GetParticleChangeF
206 isInitialised = true;
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oo
210
211 void G4DNAPTBElasticModel::ReadDiffCSFile(cons
212 cons
213 cons
214 {
215 // Method to read and save the information c
216 // This method is not yet standard.
217
218 // get the path of the G4LEDATA data folder
219 const char* path = G4FindDataDir("G4LEDATA")
220 // if it is not found then quit and print er
221 if (path == nullptr) {
222 G4Exception("G4DNAPTBElasticModel::ReadAll
223 "G4LEDATA environment variable
224 return;
225 }
226
227 // build the fullFileName path of the data f
228 std::ostringstream fullFileName;
229 fullFileName << path << "/" << file << ".dat
230
231 // open the data file
232 std::ifstream diffCrossSection(fullFileName.
233 // error if file is not there
234 std::stringstream endPath;
235 if (!diffCrossSection) {
236 endPath << "Missing data file: " << file;
237 G4Exception("G4DNAPTBElasticModel::Initial
238 endPath.str().c_str());
239 }
240
241 tValuesVec[materialName][particleName].push_
242
243 G4String line;
244
245 // read the file line by line until we reach
246 while (std::getline(diffCrossSection, line))
247 // check if the line is comment or empty
248 //
249 std::istringstream testIss(line);
250 G4String test;
251 testIss >> test;
252 // check first caracter to determine if fo
253 if (test == "#") {
254 // skip the line by beginning a new whil
255 continue;
256 }
257 // check if line is empty
258 if (line.empty()) {
259 // skip the line by beginning a new whil
260 continue;
261 }
262 //
263 // end of the check
264
265 // transform the line into a iss
266 std::istringstream iss(line);
267
268 // Variables to be filled by the input fil
269 G4double tDummy;
270 G4double eDummy;
271
272 // fill the variables with the content of
273 iss >> tDummy >> eDummy;
274
275 // SI : mandatory Vecm initialization
276
277 // Fill two vectors contained in maps of t
278 // [materialName][particleName]=vector
279 // [materialName][particleName][T]=vector
280 // to list all the incident energies (tVal
281 // file
282 //
283 // Check if we already have the current T
284 // If not then add it
285 if (tDummy != tValuesVec[materialName][par
286 // Add the current T value
287 tValuesVec[materialName][particleName].p
288 // Make it correspond to a default zero
289 eValuesVect[materialName][particleName][
290 }
291
292 // Put the differential cross section valu
293 // map
294 iss >> diffCrossSectionData[materialName][
295
296 // If the current E value (eDummy) is diff
297 // then add it to the vector
298 if (eDummy != eValuesVect[materialName][pa
299 eValuesVect[materialName][particleName][
300 }
301 }
302 }
303
304 //....oooOO0OOooo........oooOO0OOooo........oo
305
306 G4double G4DNAPTBElasticModel::CrossSectionPer
307
308
309 {
310 if (verboseLevel > 3){
311 G4cout << "Calling CrossSectionPerVolume()
312 }
313
314 // Get the name of the current particle
315 const G4String& particleName = p->GetParticl
316 const std::size_t& materialID = pMaterial->G
317
318 // set killBelowEnergy value for current mat
319 fKillBelowEnergy = fpModelData->GetLowELimit
320 // initialise the return value (cross sectio
321 G4double sigma = 0.;
322
323 // check if we are below the high energy lim
324 if (ekin < fpModelData->GetHighELimit(materi
325 // This is used to kill the particle if it
326 // If the energy is lower then we return a
327 // method will be called for sure. SampleS
328 // simulation.
329 //
330 // SI : XS must not be zero otherwise samp
331 if (ekin < fKillBelowEnergy) {
332 return DBL_MAX;
333 }
334
335 // Get the tables with the cross section d
336 auto tableData = fpModelData->GetData();
337 if ((*tableData)[materialID][p] == nullptr
338 G4Exception("G4DNAPTBElasticModel::Cross
339 "No model is registered");
340 }
341 // Retrieve the cross section value
342 sigma = (*tableData)[materialID][p]->FindV
343 }
344
345 if (verboseLevel > 2) {
346 G4cout << "_______________________________
347 G4cout << "°°° G4DNAPTBElasticModel - X
348 G4cout << "°°° Kinetic energy(eV)=" <<
349 G4cout << "°°° Cross section per molecu
350 G4cout << "°°° G4DNAPTBElasticModel - X
351 }
352
353 // Return the cross section
354 auto MolDensity =
355 (*G4DNAMolecularMaterial::Instance()->GetN
356 return sigma * MolDensity;
357 }
358
359 //....oooOO0OOooo........oooOO0OOooo........oo
360
361 void G4DNAPTBElasticModel::SampleSecondaries(s
362 c
363 c
364 G
365 {
366 if (verboseLevel > 3) {
367 G4cout << "Calling SampleSecondaries() of
368 }
369
370 G4double electronEnergy0 = aDynamicElectron-
371 const std::size_t& materialID = couple->GetM
372 auto p = aDynamicElectron->GetParticleDefini
373
374 // set killBelowEnergy value for material
375 fKillBelowEnergy = fpModelData->GetLowELimit
376
377 // If the particle (electron here) energy is
378 // simulation
379 if (electronEnergy0 < fKillBelowEnergy) {
380 fParticleChangeForGamma->SetProposedKineti
381 fParticleChangeForGamma->ProposeTrackStatu
382 fParticleChangeForGamma->ProposeLocalEnerg
383 }
384 // If we are above the kill limite and below
385 else if (electronEnergy0 >= fKillBelowEnergy
386 // Random sampling of the cosTheta
387 G4double cosTheta = fpModelData->Randomize
388
389 // Random sampling of phi
390 G4double phi = 2. * CLHEP::pi * G4UniformR
391
392 auto zVers = aDynamicElectron->GetMomentum
393 auto xVers = zVers.orthogonal();
394 auto yVers = zVers.cross(xVers);
395
396 G4double xDir = std::sqrt(1. - cosTheta *
397 G4double yDir = xDir;
398 xDir *= std::cos(phi);
399 yDir *= std::sin(phi);
400
401 // Particle direction after ModelInterface
402 G4ThreeVector zPrikeVers((xDir * xVers + y
403
404 // Give the new direction
405 fParticleChangeForGamma->ProposeMomentumDi
406
407 // Update the energy which does not change
408 fParticleChangeForGamma->SetProposedKineti
409 }
410 }
411
412 //....oooOO0OOooo........oooOO0OOooo........oo
413
414 G4double G4DNAPTBElasticModel::Theta(const G4P
415 const std
416 {
417 G4double theta = 0.;
418 G4double valueT1 = 0;
419 G4double valueT2 = 0;
420 G4double valueE21 = 0;
421 G4double valueE22 = 0;
422 G4double valueE12 = 0;
423 G4double valueE11 = 0;
424 G4double xs11 = 0;
425 G4double xs12 = 0;
426 G4double xs21 = 0;
427 G4double xs22 = 0;
428 if (p == G4Electron::ElectronDefinition()) {
429 auto t2 =
430 std::upper_bound(tValuesVec[materialID][
431 auto t1 = t2 - 1;
432
433 auto e12 = std::upper_bound(eValuesVect[ma
434 eValuesVect[ma
435 auto e11 = e12 - 1;
436
437 auto e22 = std::upper_bound(eValuesVect[ma
438 eValuesVect[ma
439 auto e21 = e22 - 1;
440
441 valueT1 = *t1;
442 valueT2 = *t2;
443 valueE21 = *e21;
444 valueE22 = *e22;
445 valueE12 = *e12;
446 valueE11 = *e11;
447
448 xs11 = diffCrossSectionData[materialID][p]
449 xs12 = diffCrossSectionData[materialID][p]
450 xs21 = diffCrossSectionData[materialID][p]
451 xs22 = diffCrossSectionData[materialID][p]
452 }
453
454 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x
455 return (0.);
456 }
457
458 theta = QuadInterpolator(valueE11, valueE12,
459 valueT2, k, integrD
460
461 return theta;
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oo
465
466 G4double G4DNAPTBElasticModel::LinLogInterpola
467
468 {
469 G4double d1 = std::log(xs1);
470 G4double d2 = std::log(xs2);
471 G4double value = std::exp(d1 + (d2 - d1) * (
472 return value;
473 }
474
475 //....oooOO0OOooo........oooOO0OOooo........oo
476
477 G4double G4DNAPTBElasticModel::LinLinInterpola
478
479 {
480 G4double d1 = xs1;
481 G4double d2 = xs2;
482 G4double value = (d1 + (d2 - d1) * (e - e1)
483 return value;
484 }
485
486 //....oooOO0OOooo........oooOO0OOooo........oo
487
488 G4double G4DNAPTBElasticModel::LogLogInterpola
489
490 {
491 G4double a = (std::log10(xs2) - std::log10(x
492 G4double b = std::log10(xs2) - a * std::log1
493 G4double sigma = a * std::log10(e) + b;
494 G4double value = (std::pow(10., sigma));
495 return value;
496 }
497
498 //....oooOO0OOooo........oooOO0OOooo........oo
499
500 G4double G4DNAPTBElasticModel::QuadInterpolato
501
502
503
504 {
505 // Log-Log
506 /*
507 G4double interpolatedvalue1 = LogLogInterpo
508 G4double interpolatedvalue2 = LogLogInterpo
509 G4double value = LogLogInterpolate(t1, t2,
510
511
512 // Lin-Log
513 G4double interpolatedvalue1 = LinLogInterpo
514 G4double interpolatedvalue2 = LinLogInterpo
515 G4double value = LinLogInterpolate(t1, t2,
516 */
517
518 // Lin-Lin
519 G4double interpolatedvalue1 = LinLinInterpol
520 G4double interpolatedvalue2 = LinLinInterpol
521 G4double value = LinLinInterpolate(t1, t2, t
522
523 return value;
524 }
525
526 //....oooOO0OOooo........oooOO0OOooo........oo
527
528 G4double G4DNAPTBElasticModel::RandomizeCosThe
529 {
530 G4double integrdiff = 0;
531 G4double uniformRand = G4UniformRand();
532 integrdiff = uniformRand;
533
534 G4double theta = 0.;
535 G4double cosTheta = 0.;
536 theta = Theta(G4Electron::ElectronDefinition
537
538 cosTheta = std::cos(theta * CLHEP::pi / 180)
539
540 return cosTheta;
541 }
542