| 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 //
27
28 #include "G4DNAChampionElasticModel.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4Exp.hh"
33
34 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35
36 using namespace std;
37
38 #define CHAMPION_VERBOSE
39
40 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
41
42 G4DNAChampionElasticModel::
43 G4DNAChampionElasticModel(const G4ParticleDefinition*, const G4String& nam) :
44 G4VEmModel(nam)
45 {
46 SetLowEnergyLimit(7.4 * eV);
47 SetHighEnergyLimit(1. * MeV);
48
49 verboseLevel = 0;
50 // Verbosity scale:
51 // 0 = nothing
52 // 1 = warning for energy non-conservation
53 // 2 = details of energy budget
54 // 3 = calculation of cross sections, file openings, sampling of atoms
55 // 4 = entering in methods
56
57 #ifdef CHAMPION_VERBOSE
58 if (verboseLevel > 0)
59 {
60 G4cout << "Champion Elastic model is constructed "
61 << G4endl
62 << "Energy range: "
63 << LowEnergyLimit() / eV << " eV - "
64 << HighEnergyLimit() / MeV << " MeV"
65 << G4endl;
66 }
67 #endif
68
69 fParticleChangeForGamma = nullptr;
70 fpMolWaterDensity = nullptr;
71 fpData = nullptr;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76 G4DNAChampionElasticModel::~G4DNAChampionElasticModel()
77 {
78 // For total cross section
79 delete fpData;
80
81 // For final state
82 eVecm.clear();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
86
87 void G4DNAChampionElasticModel::Initialise(const G4ParticleDefinition* particle,
88 const G4DataVector& /*cuts*/)
89 {
90 #ifdef CHAMPION_VERBOSE
91 if (verboseLevel > 3)
92 {
93 G4cout << "Calling G4DNAChampionElasticModel::Initialise()" << G4endl;
94 }
95 #endif
96
97 if(particle->GetParticleName() != "e-")
98 {
99 G4Exception("G4DNAChampionElasticModel::Initialise",
100 "em0002",
101 FatalException,
102 "Model not applicable to particle type.");
103 }
104
105 // Energy limits
106
107 if (LowEnergyLimit() < 7.4*eV)
108 {
109 G4cout << "G4DNAChampionElasticModel: low energy limit increased from "
110 << LowEnergyLimit()/eV << " eV to " << 7.4 << " eV"
111 << G4endl;
112 SetLowEnergyLimit(7.4*eV);
113 }
114
115 if (HighEnergyLimit() > 1.*MeV)
116 {
117 G4cout << "G4DNAChampionElasticModel: high energy limit decreased from "
118 << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
119 << G4endl;
120 SetHighEnergyLimit(1.*MeV);
121 }
122
123 if (isInitialised) { return; }
124
125 // *** ELECTRON
126 // For total cross section
127 // Reading of data files
128
129 G4double scaleFactor = 1e-16*cm*cm;
130
131 G4String fileElectron("dna/sigma_elastic_e_champion");
132
133 fpData = new G4DNACrossSectionDataSet(new G4LogLogInterpolation(),
134 eV,
135 scaleFactor );
136 fpData->LoadData(fileElectron);
137
138 // For final state
139
140 const char *path = G4FindDataDir("G4LEDATA");
141
142 if (path == nullptr)
143 {
144 G4Exception("G4ChampionElasticModel::Initialise",
145 "em0006",
146 FatalException,
147 "G4LEDATA environment variable not set.");
148 return;
149 }
150
151 std::ostringstream eFullFileName;
152 eFullFileName << path << "/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat";
153 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
154
155 if (!eDiffCrossSection)
156 {
157 G4ExceptionDescription errMsg;
158 errMsg << "Missing data file:/dna/sigmadiff_cumulated_elastic_e_champion_hp.dat; "
159 << "please use G4EMLOW6.36 and above.";
160
161 G4Exception("G4DNAChampionElasticModel::Initialise",
162 "em0003",
163 FatalException,
164 errMsg);
165 }
166
167 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
168 // Added clear for MT
169
170 eTdummyVec.clear();
171 eVecm.clear();
172 eDiffCrossSectionData.clear();
173
174 //
175
176 eTdummyVec.push_back(0.);
177
178 while(!eDiffCrossSection.eof())
179 {
180 G4double tDummy;
181 G4double eDummy;
182 eDiffCrossSection >> tDummy >> eDummy;
183
184 // SI : mandatory eVecm initialization
185
186 if (tDummy != eTdummyVec.back())
187 {
188 eTdummyVec.push_back(tDummy);
189 eVecm[tDummy].push_back(0.);
190 }
191
192 eDiffCrossSection >> eDiffCrossSectionData[tDummy][eDummy];
193
194 if (eDummy != eVecm[tDummy].back()) eVecm[tDummy].push_back(eDummy);
195 }
196
197 // End final state
198
199 #ifdef CHAMPION_VERBOSE
200 if (verboseLevel>0)
201 {
202 if (verboseLevel > 2)
203 {
204 G4cout << "Loaded cross section files for Champion Elastic model" << G4endl;
205 }
206
207 G4cout << "Champion Elastic model is initialized " << G4endl
208 << "Energy range: "
209 << LowEnergyLimit() / eV << " eV - "
210 << HighEnergyLimit() / MeV << " MeV"
211 << G4endl;
212 }
213 #endif
214
215 // Initialize water density pointer
216 G4DNAMolecularMaterial::Instance()->Initialize();
217
218 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->
219 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
220
221 fParticleChangeForGamma = GetParticleChangeForGamma();
222 isInitialised = true;
223
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227
228 G4double
229 G4DNAChampionElasticModel::
230 CrossSectionPerVolume(const G4Material* material,
231 #ifdef CHAMPION_VERBOSE
232 const G4ParticleDefinition* p,
233 #else
234 const G4ParticleDefinition*,
235 #endif
236 G4double ekin,
237 G4double,
238 G4double)
239 {
240 #ifdef CHAMPION_VERBOSE
241 if (verboseLevel > 3)
242 {
243 G4cout << "Calling CrossSectionPerVolume() of G4DNAChampionElasticModel"
244 << G4endl;
245 }
246 #endif
247
248 // Calculate total cross section for model
249
250 G4double sigma = 0.;
251 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
252
253 if (ekin <= HighEnergyLimit() && ekin >= LowEnergyLimit())
254 {
255 //SI : XS must not be zero otherwise sampling of secondaries method ignored
256 //
257 sigma = fpData->FindValue(ekin);
258 }
259
260 #ifdef CHAMPION_VERBOSE
261 if (verboseLevel > 2)
262 {
263 G4cout << "__________________________________" << G4endl;
264 G4cout << "=== G4DNAChampionElasticModel - XS INFO START" << G4endl;
265 G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : " << p->GetParticleName() << G4endl;
266 G4cout << "=== Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
267 G4cout << "=== Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
268 G4cout << "=== G4DNAChampionElasticModel - XS INFO END" << G4endl;
269 }
270 #endif
271
272 return sigma*waterDensity;
273 }
274
275 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
276
277 void G4DNAChampionElasticModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
278 const G4MaterialCutsCouple* /*couple*/,
279 const G4DynamicParticle* aDynamicElectron,
280 G4double,
281 G4double)
282 {
283
284 #ifdef CHAMPION_VERBOSE
285 if (verboseLevel > 3)
286 {
287 G4cout << "Calling SampleSecondaries() of G4DNAChampionElasticModel" << G4endl;
288 }
289 #endif
290
291 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
292
293 G4double cosTheta = RandomizeCosTheta(electronEnergy0);
294
295 G4double phi = 2. * pi * G4UniformRand();
296
297 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
298 G4ThreeVector xVers = zVers.orthogonal();
299 G4ThreeVector yVers = zVers.cross(xVers);
300
301 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
302 G4double yDir = xDir;
303 xDir *= std::cos(phi);
304 yDir *= std::sin(phi);
305
306 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
307
308 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
309
310 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
311
312 }
313
314 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316 G4double G4DNAChampionElasticModel::Theta(G4double k,
317 G4double integrDiff)
318 {
319 G4double theta = 0.;
320 G4double valueT1 = 0;
321 G4double valueT2 = 0;
322 G4double valueE21 = 0;
323 G4double valueE22 = 0;
324 G4double valueE12 = 0;
325 G4double valueE11 = 0;
326 G4double xs11 = 0;
327 G4double xs12 = 0;
328 G4double xs21 = 0;
329 G4double xs22 = 0;
330
331 auto t2 = std::upper_bound(eTdummyVec.begin(),
332 eTdummyVec.end(), k);
333 auto t1 = t2 - 1;
334
335 auto e12 = std::upper_bound(eVecm[(*t1)].begin(),
336 eVecm[(*t1)].end(),
337 integrDiff);
338 auto e11 = e12 - 1;
339
340 auto e22 = std::upper_bound(eVecm[(*t2)].begin(),
341 eVecm[(*t2)].end(),
342 integrDiff);
343 auto e21 = e22 - 1;
344
345 valueT1 = *t1;
346 valueT2 = *t2;
347 valueE21 = *e21;
348 valueE22 = *e22;
349 valueE12 = *e12;
350 valueE11 = *e11;
351
352 xs11 = eDiffCrossSectionData[valueT1][valueE11];
353 xs12 = eDiffCrossSectionData[valueT1][valueE12];
354 xs21 = eDiffCrossSectionData[valueT2][valueE21];
355 xs22 = eDiffCrossSectionData[valueT2][valueE22];
356
357 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && xs22 == 0) return (0.);
358
359 theta = QuadInterpolator(valueE11, valueE12, valueE21, valueE22, xs11, xs12,
360 xs21, xs22, valueT1, valueT2, k, integrDiff);
361
362 return theta;
363 }
364
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366
367 G4double G4DNAChampionElasticModel::LinLogInterpolate(G4double e1,
368 G4double e2,
369 G4double e,
370 G4double xs1,
371 G4double xs2)
372 {
373 G4double d1 = std::log(xs1);
374 G4double d2 = std::log(xs2);
375 G4double value = G4Exp(d1 + (d2 - d1) * (e - e1) / (e2 - e1));
376 return value;
377 }
378
379 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
380
381 G4double G4DNAChampionElasticModel::LinLinInterpolate(G4double e1,
382 G4double e2,
383 G4double e,
384 G4double xs1,
385 G4double xs2)
386 {
387 G4double d1 = xs1;
388 G4double d2 = xs2;
389 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
390 return value;
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
394
395 G4double G4DNAChampionElasticModel::LogLogInterpolate(G4double e1,
396 G4double e2,
397 G4double e,
398 G4double xs1,
399 G4double xs2)
400 {
401 G4double a = (std::log10(xs2) - std::log10(xs1))
402 / (std::log10(e2) - std::log10(e1));
403 G4double b = std::log10(xs2) - a * std::log10(e2);
404 G4double sigma = a * std::log10(e) + b;
405 G4double value = (std::pow(10., sigma));
406 return value;
407 }
408
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
410
411 G4double G4DNAChampionElasticModel::QuadInterpolator(G4double e11,
412 G4double e12,
413 G4double e21,
414 G4double e22,
415 G4double xs11,
416 G4double xs12,
417 G4double xs21,
418 G4double xs22,
419 G4double t1,
420 G4double t2,
421 G4double t,
422 G4double e)
423 {
424 // Log-Log
425 /*
426 G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
427 G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
428 G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
429
430
431 // Lin-Log
432 G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
433 G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
434 G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
435 */
436
437 // Lin-Lin
438 G4double interpolatedvalue1 = LinLinInterpolate(e11, e12, e, xs11, xs12);
439 G4double interpolatedvalue2 = LinLinInterpolate(e21, e22, e, xs21, xs22);
440 G4double value = LinLinInterpolate(t1, t2, t, interpolatedvalue1,
441 interpolatedvalue2);
442
443 return value;
444 }
445
446 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
447
448 G4double G4DNAChampionElasticModel::RandomizeCosTheta(G4double k)
449 {
450
451 G4double integrdiff = 0;
452 G4double uniformRand = G4UniformRand();
453 integrdiff = uniformRand;
454
455 G4double theta = 0.;
456 G4double cosTheta = 0.;
457 theta = Theta(k / eV, integrdiff);
458
459 cosTheta = std::cos(theta * pi / 180);
460
461 return cosTheta;
462 }
463
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465
466 void G4DNAChampionElasticModel::SetKillBelowThreshold(G4double)
467 {
468 G4ExceptionDescription errMsg;
469 errMsg << "The method G4DNAChampionElasticModel::SetKillBelowThreshold is deprecated";
470
471 G4Exception("G4DNAChampionElasticModel::SetKillBelowThreshold",
472 "deprecated",
473 JustWarning,
474 errMsg);
475 }
476