| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Models come from
28 // M. Bug et al, Rad. Phys and Chem. 130, 459-479 (2017)
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(const G4String& applyToMaterial,
40 const G4ParticleDefinition*, const G4String& nam)
41 : G4VDNAModel(nam, applyToMaterial)
42 {
43 if (verboseLevel > 0) {
44 G4cout << "PTB Elastic model is constructed : " << G4endl;
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_WATER", false);
51 fpBackbone_THF = G4Material::GetMaterial("backbone_THF", false);
52 fpCytosine_PY = G4Material::GetMaterial("cytosine_PY", false);
53 fpThymine_PY = G4Material::GetMaterial("thymine_PY", false);
54 fpAdenine_PU = G4Material::GetMaterial("adenine_PU", false);
55 fpBackbone_TMP = G4Material::GetMaterial("backbone_TMP", false);
56 fpGuanine_PU = G4Material::GetMaterial("guanine_PU", false);
57 fpN2 = G4Material::GetMaterial("N2", false);
58 }
59
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
61
62 void G4DNAPTBElasticModel::Initialise(const G4ParticleDefinition* particle,
63 const G4DataVector& /*cuts*/)
64 {
65 if (isInitialised) {
66 return;
67 }
68 if (verboseLevel > 3)
69 {
70 G4cout << "Calling G4DNAPTBElasticModel::Initialise()" << G4endl;
71 }
72 if (particle != G4Electron::ElectronDefinition()) {
73 std::ostringstream oss;
74 oss << " Model is not applied for this particle " << particle->GetParticleName();
75 G4Exception("G4DNAPTBElasticModel::G4DNAPTBElasticModel", "PTB001", FatalException,
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/sigma_elastic_e-_PTB_N2",
88 "dna/sigmadiff_cumulated_elastic_e-_PTB_N2", scaleFactor);
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/sigma_elastic_e-_PTB_THF",
97 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF", scaleFactor);
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/sigma_elastic_e-_PTB_PY",
105 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY", scaleFactor);
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/sigma_elastic_e-_PTB_PU",
113 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU", scaleFactor);
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/sigma_elastic_e-_PTB_TMP",
121 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP", scaleFactor);
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/sigma_elastic_e_champion",
129 "dna/sigmadiff_cumulated_elastic_e_champion", scaleFactor);
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/sigma_elastic_e-_PTB_THF",
138 "dna/sigmadiff_cumulated_elastic_e-_PTB_THF", scaleFactor * 33. / 30);
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/sigma_elastic_e-_PTB_PY",
146 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY", scaleFactor * 42. / 30);
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/sigma_elastic_e-_PTB_PY",
154 "dna/sigmadiff_cumulated_elastic_e-_PTB_PY", scaleFactor * 48. / 30);
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/sigma_elastic_e-_PTB_PU",
162 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU", scaleFactor * 50. / 44);
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/sigma_elastic_e-_PTB_PU",
169 "dna/sigmadiff_cumulated_elastic_e-_PTB_PU", scaleFactor * 56. / 44);
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/sigma_elastic_e-_PTB_TMP",
177 "dna/sigmadiff_cumulated_elastic_e-_PTB_TMP", scaleFactor * 33. / 50);
178 SetLowELimit(index, particle, 10 * eV);
179 SetHighELimit(index, particle, 1 * keV);
180 }
181
182 if (!G4DNAMaterialManager::Instance()->IsLocked()) {
183 // Load the data
184 LoadCrossSectionData(particle);
185 G4DNAMaterialManager::Instance()->SetMasterDataModel(DNAModelType::fDNAElastics, this);
186 fpModelData = this;
187 }
188 else {
189 auto dataModel = dynamic_cast<G4DNAPTBElasticModel*>(
190 G4DNAMaterialManager::Instance()->GetModel(DNAModelType::fDNAElastics));
191 if (dataModel == nullptr) {
192 G4cout << "G4DNAPTBElasticModel::Initialise:: not good modelData" << G4endl;
193 G4Exception("G4DNAPTBElasticModel::Initialise", "PTB0006", FatalException,
194 "not good modelData");
195 }
196 else {
197 fpModelData = dataModel;
198 }
199 }
200
201 if (verboseLevel > 2) {
202 G4cout << "Loaded cross section files for PTB Elastic model" << G4endl;
203 }
204
205 fParticleChangeForGamma = GetParticleChangeForGamma();
206 isInitialised = true;
207 }
208
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210
211 void G4DNAPTBElasticModel::ReadDiffCSFile(const std::size_t& materialName,
212 const G4ParticleDefinition* particleName,
213 const G4String& file, const G4double&)
214 {
215 // Method to read and save the information contained within the differential cross section files.
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 error message
221 if (path == nullptr) {
222 G4Exception("G4DNAPTBElasticModel::ReadAllDiffCSFiles", "em0006", FatalException,
223 "G4LEDATA environment variable not set.");
224 return;
225 }
226
227 // build the fullFileName path of the data file
228 std::ostringstream fullFileName;
229 fullFileName << path << "/" << file << ".dat";
230
231 // open the data file
232 std::ifstream diffCrossSection(fullFileName.str().c_str());
233 // error if file is not there
234 std::stringstream endPath;
235 if (!diffCrossSection) {
236 endPath << "Missing data file: " << file;
237 G4Exception("G4DNAPTBElasticModel::Initialise", "em0003", FatalException,
238 endPath.str().c_str());
239 }
240
241 tValuesVec[materialName][particleName].push_back(0.);
242
243 G4String line;
244
245 // read the file line by line until we reach the end of file point
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 following information is data or comments
253 if (test == "#") {
254 // skip the line by beginning a new while loop.
255 continue;
256 }
257 // check if line is empty
258 if (line.empty()) {
259 // skip the line by beginning a new while loop.
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 file
269 G4double tDummy;
270 G4double eDummy;
271
272 // fill the variables with the content of the line
273 iss >> tDummy >> eDummy;
274
275 // SI : mandatory Vecm initialization
276
277 // Fill two vectors contained in maps of types:
278 // [materialName][particleName]=vector
279 // [materialName][particleName][T]=vector
280 // to list all the incident energies (tValues) and all the output energies (eValues) within the
281 // file
282 //
283 // Check if we already have the current T value in the vector.
284 // If not then add it
285 if (tDummy != tValuesVec[materialName][particleName].back()) {
286 // Add the current T value
287 tValuesVec[materialName][particleName].push_back(tDummy);
288 // Make it correspond to a default zero E value
289 eValuesVect[materialName][particleName][tDummy].push_back(0.);
290 }
291
292 // Put the differential cross section value of the input file within the diffCrossSectionData
293 // map
294 iss >> diffCrossSectionData[materialName][particleName][tDummy][eDummy];
295
296 // If the current E value (eDummy) is different from the one already registered in the eVector
297 // then add it to the vector
298 if (eDummy != eValuesVect[materialName][particleName][tDummy].back()) {
299 eValuesVect[materialName][particleName][tDummy].push_back(eDummy);
300 }
301 }
302 }
303
304 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305
306 G4double G4DNAPTBElasticModel::CrossSectionPerVolume(const G4Material* pMaterial,
307 const G4ParticleDefinition* p, G4double ekin,
308 G4double /*emin*/, G4double /*emax*/)
309 {
310 if (verboseLevel > 3){
311 G4cout << "Calling CrossSectionPerVolume() of G4DNAPTBElasticModel" << G4endl;
312 }
313
314 // Get the name of the current particle
315 const G4String& particleName = p->GetParticleName();
316 const std::size_t& materialID = pMaterial->GetIndex();
317
318 // set killBelowEnergy value for current material
319 fKillBelowEnergy = fpModelData->GetLowELimit(materialID, p);
320 // initialise the return value (cross section) to zero
321 G4double sigma = 0.;
322
323 // check if we are below the high energy limit
324 if (ekin < fpModelData->GetHighELimit(materialID, p)) {
325 // This is used to kill the particle if its kinetic energy is below fKillBelowEnergy.
326 // If the energy is lower then we return a maximum cross section and thus the SampleSecondaries
327 // method will be called for sure. SampleSecondaries will remove the particle from the
328 // simulation.
329 //
330 // SI : XS must not be zero otherwise sampling of secondaries method ignored
331 if (ekin < fKillBelowEnergy) {
332 return DBL_MAX;
333 }
334
335 // Get the tables with the cross section data
336 auto tableData = fpModelData->GetData();
337 if ((*tableData)[materialID][p] == nullptr) {
338 G4Exception("G4DNAPTBElasticModel::CrossSectionPerVolume", "em00236", FatalException,
339 "No model is registered");
340 }
341 // Retrieve the cross section value
342 sigma = (*tableData)[materialID][p]->FindValue(ekin);
343 }
344
345 if (verboseLevel > 2) {
346 G4cout << "__________________________________" << G4endl;
347 G4cout << "°°° G4DNAPTBElasticModel - XS INFO START" << G4endl;
348 G4cout << "°°° Kinetic energy(eV)=" << ekin / eV << " particle : " << particleName << G4endl;
349 G4cout << "°°° Cross section per molecule (cm^2)=" << sigma / cm / cm << G4endl;
350 G4cout << "°°° G4DNAPTBElasticModel - XS INFO END" << G4endl;
351 }
352
353 // Return the cross section
354 auto MolDensity =
355 (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(pMaterial))[materialID];
356 return sigma * MolDensity;
357 }
358
359 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
360
361 void G4DNAPTBElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
362 const G4MaterialCutsCouple* couple,
363 const G4DynamicParticle* aDynamicElectron,
364 G4double /*tmin*/, G4double /*tmax*/)
365 {
366 if (verboseLevel > 3) {
367 G4cout << "Calling SampleSecondaries() of G4DNAPTBElasticModel" << G4endl;
368 }
369
370 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
371 const std::size_t& materialID = couple->GetMaterial()->GetIndex();
372 auto p = aDynamicElectron->GetParticleDefinition();
373
374 // set killBelowEnergy value for material
375 fKillBelowEnergy = fpModelData->GetLowELimit(materialID, p);
376
377 // If the particle (electron here) energy is below the kill limit then we remove it from the
378 // simulation
379 if (electronEnergy0 < fKillBelowEnergy) {
380 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
381 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
382 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
383 }
384 // If we are above the kill limite and below the high limit then we proceed
385 else if (electronEnergy0 >= fKillBelowEnergy && electronEnergy0 < GetHighELimit(materialID, p)) {
386 // Random sampling of the cosTheta
387 G4double cosTheta = fpModelData->RandomizeCosTheta(electronEnergy0, materialID);
388
389 // Random sampling of phi
390 G4double phi = 2. * CLHEP::pi * G4UniformRand();
391
392 auto zVers = aDynamicElectron->GetMomentumDirection();
393 auto xVers = zVers.orthogonal();
394 auto yVers = zVers.cross(xVers);
395
396 G4double xDir = std::sqrt(1. - cosTheta * 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 + yDir * yVers + cosTheta * zVers));
403
404 // Give the new direction
405 fParticleChangeForGamma->ProposeMomentumDirection(zPrikeVers.unit());
406
407 // Update the energy which does not change here
408 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
409 }
410 }
411
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413
414 G4double G4DNAPTBElasticModel::Theta(const G4ParticleDefinition* p, G4double k, G4double integrDiff,
415 const std::size_t& materialID)
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][p].begin(), tValuesVec[materialID][p].end(), k);
431 auto t1 = t2 - 1;
432
433 auto e12 = std::upper_bound(eValuesVect[materialID][p][(*t1)].begin(),
434 eValuesVect[materialID][p][(*t1)].end(), integrDiff);
435 auto e11 = e12 - 1;
436
437 auto e22 = std::upper_bound(eValuesVect[materialID][p][(*t2)].begin(),
438 eValuesVect[materialID][p][(*t2)].end(), integrDiff);
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][valueT1][valueE11];
449 xs12 = diffCrossSectionData[materialID][p][valueT1][valueE12];
450 xs21 = diffCrossSectionData[materialID][p][valueT2][valueE21];
451 xs22 = diffCrossSectionData[materialID][p][valueT2][valueE22];
452 }
453
454 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) {
455 return (0.);
456 }
457
458 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, xs21, xs22, valueT1,
459 valueT2, k, integrDiff);
460
461 return theta;
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465
466 G4double G4DNAPTBElasticModel::LinLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1,
467 G4double xs2)
468 {
469 G4double d1 = std::log(xs1);
470 G4double d2 = std::log(xs2);
471 G4double value = std::exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
472 return value;
473 }
474
475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
476
477 G4double G4DNAPTBElasticModel::LinLinInterpolate(G4double e1, G4double e2, G4double e, G4double xs1,
478 G4double xs2)
479 {
480 G4double d1 = xs1;
481 G4double d2 = xs2;
482 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
483 return value;
484 }
485
486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487
488 G4double G4DNAPTBElasticModel::LogLogInterpolate(G4double e1, G4double e2, G4double e, G4double xs1,
489 G4double xs2)
490 {
491 G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
492 G4double b = std::log10(xs2) - a * std::log10(e2);
493 G4double sigma = a * std::log10(e) + b;
494 G4double value = (std::pow(10., sigma));
495 return value;
496 }
497
498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
499
500 G4double G4DNAPTBElasticModel::QuadInterpolator(G4double e11, G4double e12, G4double e21,
501 G4double e22, G4double xs11, G4double xs12,
502 G4double xs21, G4double xs22, G4double t1,
503 G4double t2, G4double t, G4double e)
504 {
505 // Log-Log
506 /*
507 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
508 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
509 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
510
511
512 // Lin-Log
513 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
514 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
515 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
516 */
517
518 // Lin-Lin
519 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
520 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
521 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
522
523 return value;
524 }
525
526 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
527
528 G4double G4DNAPTBElasticModel::RandomizeCosTheta(const G4double& k, const std::size_t& materialID)
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(), k / eV, integrdiff, materialID);
537
538 cosTheta = std::cos(theta * CLHEP::pi / 180);
539
540 return cosTheta;
541 }
542