| 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 // CPA100 elastic model class for electrons
27 //
28 // Based on the work of M. Terrissol and M. C. Bordage
29 //
30 // Users are requested to cite the following papers:
31 // - M. Terrissol, A. Baudre, Radiat. Prot. Dosim. 31 (1990) 175-177
32 // - M.C. Bordage, J. Bordes, S. Edel, M. Terrissol, X. Franceries,
33 // M. Bardies, N. Lampe, S. Incerti, Phys. Med. 32 (2016) 1833-1840
34 //
35 // Authors of this class:
36 // M.C. Bordage, M. Terrissol, S. Edel, J. Bordes, S. Incerti
37 //
38 // 15.01.2014: creation
39 //
40 // Based on the study by S. Zein et. al. Nucl. Inst. Meth. B 488 (2021) 70-82
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::G4DNACPA100ElasticModel(const G4ParticleDefinition*, const G4String& nam)
51 : G4VDNAModel(nam, "all")
52 {
53 fpGuanine = G4Material::GetMaterial("G4_GUANINE", false);
54 fpG4_WATER = G4Material::GetMaterial("G4_WATER", false);
55 fpDeoxyribose = G4Material::GetMaterial("G4_DEOXYRIBOSE", false);
56 fpCytosine = G4Material::GetMaterial("G4_CYTOSINE", false);
57 fpThymine = G4Material::GetMaterial("G4_THYMINE", false);
58 fpAdenine = G4Material::GetMaterial("G4_ADENINE", false);
59 fpPhosphate = G4Material::GetMaterial("G4_PHOSPHORIC_ACID", false);
60 fpParticle = G4Electron::ElectronDefinition();
61 }
62
63 void G4DNACPA100ElasticModel::Initialise(const G4ParticleDefinition* p,
64 const G4DataVector& /*cuts*/)
65 {
66 if (isInitialised) {
67 return;
68 }
69 if (verboseLevel > 3) {
70 G4cout << "Calling G4DNACPA100ExcitationModel::Initialise()" << G4endl;
71 }
72
73 if (!G4DNAMaterialManager::Instance()->IsLocked()) {
74 if (p != fpParticle) {
75 std::ostringstream oss;
76 oss << " Model is not applied for this particle " << p->GetParticleName();
77 G4Exception("G4DNACPA100ElasticModel::G4DNACPA100ElasticModel", "CPA001", FatalException,
78 oss.str().c_str());
79 }
80
81 const char* path = G4FindDataDir("G4LEDATA");
82
83 if (path == nullptr) {
84 G4Exception("G4DNACPA100ElasticModel::Initialise", "em0006", FatalException,
85 "G4LEDATA environment variable not set.");
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_elastic_e_cpa100",
94 "dna/sigmadiff_cumulated_elastic_e_cpa100", 1e-20 * m * m);
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_elastic_e_cpa100_guanine",
102 "dna/sigmadiff_cumulated_elastic_e_cpa100_guanine", 1 * cm * cm);
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_elastic_e_cpa100_deoxyribose",
110 "dna/sigmadiff_cumulated_elastic_e_cpa100_deoxyribose", 1 * cm * cm);
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_elastic_e_cpa100_cytosine",
118 "dna/sigmadiff_cumulated_elastic_e_cpa100_cytosine", 1 * cm * cm);
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_elastic_e_cpa100_thymine",
126 "dna/sigmadiff_cumulated_elastic_e_cpa100_thymine", 1 * cm * cm);
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_elastic_e_cpa100_adenine",
134 "dna/sigmadiff_cumulated_elastic_e_cpa100_adenine", 1 * cm * cm);
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_elastic_e_cpa100_phosphoric_acid",
142 "dna/sigmadiff_cumulated_elastic_e_cpa100_phosphoric_acid", 1 * cm * cm);
143 SetLowELimit(index, p, 11 * eV);
144 SetHighELimit(index, p, 1 * MeV);
145 }
146
147 // Load data
148 LoadCrossSectionData(p);
149 G4DNAMaterialManager::Instance()->SetMasterDataModel(DNAModelType::fDNAElastics, this);
150 fpModelData = this;
151 }
152 else {
153 auto dataModel = dynamic_cast<G4DNACPA100ElasticModel*>(
154 G4DNAMaterialManager::Instance()->GetModel(DNAModelType::fDNAElastics));
155 if (dataModel == nullptr) {
156 G4cout << "G4DNACPA100ElasticModel::CrossSectionPerVolume:: not good modelData" << G4endl;
157 G4Exception("G4DNACPA100ElasticModel::CrossSectionPerVolume", "em004", FatalException,
158 "no modelData is registered");
159 }
160 else {
161 fpModelData = dataModel;
162 }
163 }
164
165 fParticleChangeForGamma = GetParticleChangeForGamma();
166 isInitialised = true;
167 }
168
169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
170
171 G4double G4DNACPA100ElasticModel::CrossSectionPerVolume(const G4Material* pMaterial,
172 const G4ParticleDefinition* p,
173 G4double ekin, G4double, G4double)
174 {
175 // Get the name of the current particle
176 const G4String& particleName = p->GetParticleName();
177 auto materialID = pMaterial->GetIndex();
178
179 // set killBelowEnergy value for current material
180 fKillBelowEnergy = fpModelData->GetLowELimit(materialID, p);
181
182 G4double sigma = 0.;
183
184 if (ekin < fpModelData->GetHighELimit(materialID, p)) {
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::CrossSectionPerVolume", "em00236", FatalException,
193 "No model is registered");
194 }
195 sigma = (*tableData)[materialID][p]->FindValue(ekin);
196 }
197
198 if (verboseLevel > 2) {
199 auto MolDensity =
200 (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(pMaterial))[materialID];
201 G4cout << "__________________________________" << G4endl;
202 G4cout << "°°° G4DNACPA100ElasticModel - XS INFO START" << G4endl;
203 G4cout << "°°° Kinetic energy(eV)=" << ekin / eV << " particle : " << particleName << G4endl;
204 G4cout << "°°° lowLim (eV) = " << GetLowELimit(materialID, p) / eV
205 << " highLim (eV) : " << GetHighELimit(materialID, p) / eV << G4endl;
206 G4cout << "°°° Materials = " << (*G4Material::GetMaterialTable())[materialID]->GetName()
207 << G4endl;
208 G4cout << "°°° Cross section per molecule (cm^2)=" << sigma / cm / cm << G4endl;
209 G4cout << "°°° Cross section per Phosphate molecule (cm^-1)=" << sigma * MolDensity / (1. / cm)
210 << G4endl;
211 G4cout << "°°° G4DNACPA100ElasticModel - XS INFO END" << G4endl;
212 }
213
214 auto MolDensity =
215 (*G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(pMaterial))[materialID];
216 return sigma * MolDensity;
217 }
218
219 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220
221 void G4DNACPA100ElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
222 const G4MaterialCutsCouple* couple,
223 const G4DynamicParticle* aDynamicElectron, G4double,
224 G4double)
225 {
226 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
227 auto materialID = couple->GetMaterial()->GetIndex();
228 auto p = aDynamicElectron->GetParticleDefinition();
229
230 if (p != fpParticle) {
231 G4Exception("G4DNACPA100ElasticModel::SampleSecondaries", "em00436", FatalException,
232 "This particle is not applied for this model");
233 }
234 if (electronEnergy0 < fKillBelowEnergy) {
235 return;
236 }
237 G4double cosTheta = fpModelData->RandomizeCosTheta(electronEnergy0, materialID);
238 G4double phi = 2. * CLHEP::pi * G4UniformRand();
239
240 const G4ThreeVector& zVers = aDynamicElectron->GetMomentumDirection();
241
242 G4double CT1, ST1, CF1, SF1, CT2, ST2, CF2, SF2;
243 G4double sinTheta = std::sqrt(1 - cosTheta * 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 * G4UniformRand());
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 * SF2, CT2);
272
273 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
274
275 auto EnergyDeposit = fpModelData->GetElasticLevel(materialID) * (1. - cosTheta) * electronEnergy0;
276 fParticleChangeForGamma->ProposeLocalEnergyDeposit(EnergyDeposit);
277 if (statCode) {
278 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
279 }
280 else {
281 auto newEnergy = electronEnergy0 - EnergyDeposit;
282 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
283 }
284 }
285
286 G4double G4DNACPA100ElasticModel::Theta(const G4ParticleDefinition* p, G4double k,
287 G4double integrDiff, const std::size_t& materialID)
288 {
289 G4double theta, valueT1, valueT2, valueE21, valueE22, valueE12, valueE11;
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][p].begin(), tValuesVec[materialID][p].end(), k);
300 auto t1 = t2 - 1;
301
302 auto e12 = std::upper_bound(eValuesVect[materialID][p][(*t1)].begin(),
303 eValuesVect[materialID][p][(*t1)].end(), integrDiff);
304 auto e11 = e12 - 1;
305
306 auto e22 = std::upper_bound(eValuesVect[materialID][p][(*t2)].begin(),
307 eValuesVect[materialID][p][(*t2)].end(), integrDiff);
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][valueT1][valueE11];
318 xs12 = diffCrossSectionData[materialID][p][valueT1][valueE12];
319 xs21 = diffCrossSectionData[materialID][p][valueT2][valueE21];
320 xs22 = diffCrossSectionData[materialID][p][valueT2][valueE22];
321 }
322
323 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) {
324 return (0.);
325 }
326
327 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12, xs21, xs22, valueT1,
328 valueT2, k, integrDiff);
329
330 return theta;
331 }
332
333 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
334
335 G4double G4DNACPA100ElasticModel::LinLogInterpolate(G4double e1, G4double e2, G4double e,
336 G4double xs1, G4double xs2)
337 {
338 G4double d1 = std::log(xs1);
339 G4double d2 = std::log(xs2);
340 G4double value = std::exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
341 return value;
342 }
343
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345
346 G4double G4DNACPA100ElasticModel::LinLinInterpolate(G4double e1, G4double e2, G4double e,
347 G4double xs1, G4double xs2)
348 {
349 G4double d1 = xs1;
350 G4double d2 = xs2;
351 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
352 return value;
353 }
354
355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
356
357 G4double G4DNACPA100ElasticModel::LogLogInterpolate(G4double e1, G4double e2, G4double e,
358 G4double xs1, G4double xs2)
359 {
360 G4double a = (std::log10(xs2) - std::log10(xs1)) / (std::log10(e2) - std::log10(e1));
361 G4double b = std::log10(xs2) - a * std::log10(e2);
362 G4double sigma = a * std::log10(e) + b;
363 G4double value = (std::pow(10., sigma));
364 return value;
365 }
366
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368
369 G4double G4DNACPA100ElasticModel::QuadInterpolator(G4double e11, G4double e12, G4double e21,
370 G4double e22, G4double xs11, G4double xs12,
371 G4double xs21, G4double xs22, G4double t1,
372 G4double t2, G4double t, G4double e)
373 {
374 // Log-Log
375 /*
376 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
377 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
378 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
379
380
381 // Lin-Log
382 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
383 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
384 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
385 */
386
387 // Lin-Lin
388 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
389 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
390 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
391
392 return value;
393 }
394
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
396
397 G4double G4DNACPA100ElasticModel::RandomizeCosTheta(G4double k, const std::size_t& materialID)
398 {
399 G4double integrdiff = 0; // PROBABILITY between 0 and 1.
400 G4double uniformRand = G4UniformRand();
401 integrdiff = uniformRand;
402 G4double cosTheta = 0.;
403 cosTheta = 1 - Theta(G4Electron::ElectronDefinition(), k / eV, integrdiff, materialID);
404 // cosTheta = std::cos(theta * CLHEP::pi / 180); ???
405 return cosTheta;
406 }
407 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
408
409 void G4DNACPA100ElasticModel::ReadDiffCSFile(const std::size_t& materialName,
410 const G4ParticleDefinition* particleName,
411 const G4String& file, const G4double&)
412 {
413 const char* path = G4FindDataDir("G4LEDATA");
414 if (path == nullptr) {
415 G4Exception("G4DNACPA100ElasticModel::ReadAllDiffCSFiles", "em0006", FatalException,
416 "G4LEDATA environment variable not set.");
417 return;
418 }
419
420 std::ostringstream fullFileName;
421 fullFileName << path << "/" << file << ".dat";
422
423 std::ifstream diffCrossSection(fullFileName.str().c_str());
424 // error if file is not there
425 std::stringstream endPath;
426 if (!diffCrossSection) {
427 endPath << "Missing data file: " << file;
428 G4Exception("G4DNACPA100ElasticModel::Initialise", "em0003", FatalException,
429 endPath.str().c_str());
430 }
431
432 tValuesVec[materialName][particleName].push_back(0.);
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][particleName].back()) {
455 // Add the current T value
456 tValuesVec[materialName][particleName].push_back(tDummy);
457 // Make it correspond to a default zero E value
458 eValuesVect[materialName][particleName][tDummy].push_back(0.);
459 }
460 iss >> diffCrossSectionData[materialName][particleName][tDummy][eDummy];
461
462 if (eDummy != eValuesVect[materialName][particleName][tDummy].back()) {
463 eValuesVect[materialName][particleName][tDummy].push_back(eDummy);
464 }
465 }
466 }
467