| 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 //
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........oo
35
36 using namespace std;
37
38 #define CHAMPION_VERBOSE
39
40 //....oooOO0OOooo........oooOO0OOooo........oo
41
42 G4DNAChampionElasticModel::
43 G4DNAChampionElasticModel(const G4ParticleDefi
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 o
55 // 4 = entering in methods
56
57 #ifdef CHAMPION_VERBOSE
58 if (verboseLevel > 0)
59 {
60 G4cout << "Champion Elastic model is const
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........oo
75
76 G4DNAChampionElasticModel::~G4DNAChampionElast
77 {
78 // For total cross section
79 delete fpData;
80
81 // For final state
82 eVecm.clear();
83 }
84
85 //....oooOO0OOooo........oooOO0OOooo........oo
86
87 void G4DNAChampionElasticModel::Initialise(con
88 con
89 {
90 #ifdef CHAMPION_VERBOSE
91 if (verboseLevel > 3)
92 {
93 G4cout << "Calling G4DNAChampionElasticMod
94 }
95 #endif
96
97 if(particle->GetParticleName() != "e-")
98 {
99 G4Exception("G4DNAChampionElasticModel::In
100 "em0002",
101 FatalException,
102 "Model not applicable to parti
103 }
104
105 // Energy limits
106
107 if (LowEnergyLimit() < 7.4*eV)
108 {
109 G4cout << "G4DNAChampionElasticModel: low
110 << LowEnergyLimit()/eV << " eV to "
111 << G4endl;
112 SetLowEnergyLimit(7.4*eV);
113 }
114
115 if (HighEnergyLimit() > 1.*MeV)
116 {
117 G4cout << "G4DNAChampionElasticModel: high
118 << HighEnergyLimit()/MeV << " MeV t
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_c
132
133 fpData = new G4DNACrossSectionDataSet(new G4
134 eV,
135 scaleF
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::Initi
145 "em0006",
146 FatalException,
147 "G4LEDATA environment variable
148 return;
149 }
150
151 std::ostringstream eFullFileName;
152 eFullFileName << path << "/dna/sigmadiff_cum
153 std::ifstream eDiffCrossSection(eFullFileNam
154
155 if (!eDiffCrossSection)
156 {
157 G4ExceptionDescription errMsg;
158 errMsg << "Missing data file:/dna/sigmadif
159 << "please use G4EMLOW6.36 and abov
160
161 G4Exception("G4DNAChampionElasticModel::In
162 "em0003",
163 FatalException,
164 errMsg);
165 }
166
167 // March 25th, 2014 - Vaclav Stepan, Sebasti
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
193
194 if (eDummy != eVecm[tDummy].back()) eVecm[
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 fo
205 }
206
207 G4cout << "Champion Elastic model is initi
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()->Initiali
217
218 fpMolWaterDensity = G4DNAMolecularMaterial::
219 GetNumMolPerVolTableFor(G4Material::GetMat
220
221 fParticleChangeForGamma = GetParticleChangeF
222 isInitialised = true;
223
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oo
227
228 G4double
229 G4DNAChampionElasticModel::
230 CrossSectionPerVolume(const G4Material* materi
231 #ifdef CHAMPION_VERBOSE
232 const G4ParticleDefiniti
233 #else
234 const G4ParticleDefiniti
235 #endif
236 G4double ekin,
237 G4double,
238 G4double)
239 {
240 #ifdef CHAMPION_VERBOSE
241 if (verboseLevel > 3)
242 {
243 G4cout << "Calling CrossSectionPerVolume()
244 << G4endl;
245 }
246 #endif
247
248 // Calculate total cross section for model
249
250 G4double sigma = 0.;
251 G4double waterDensity = (*fpMolWaterDensity)
252
253 if (ekin <= HighEnergyLimit() && ekin >= Low
254 {
255 //SI : XS must not be zero otherwise sam
256 //
257 sigma = fpData->FindValue(ekin);
258 }
259
260 #ifdef CHAMPION_VERBOSE
261 if (verboseLevel > 2)
262 {
263 G4cout << "_______________________________
264 G4cout << "=== G4DNAChampionElasticModel -
265 G4cout << "=== Kinetic energy(eV)=" << eki
266 G4cout << "=== Cross section per water mol
267 G4cout << "=== Cross section per water mol
268 G4cout << "=== G4DNAChampionElasticModel -
269 }
270 #endif
271
272 return sigma*waterDensity;
273 }
274
275 //....oooOO0OOooo........oooOO0OOooo........oo
276
277 void G4DNAChampionElasticModel::SampleSecondar
278
279
280
281
282 {
283
284 #ifdef CHAMPION_VERBOSE
285 if (verboseLevel > 3)
286 {
287 G4cout << "Calling SampleSecondaries() of
288 }
289 #endif
290
291 G4double electronEnergy0 = aDynamicElectron-
292
293 G4double cosTheta = RandomizeCosTheta(electr
294
295 G4double phi = 2. * pi * G4UniformRand();
296
297 G4ThreeVector zVers = aDynamicElectron->GetM
298 G4ThreeVector xVers = zVers.orthogonal();
299 G4ThreeVector yVers = zVers.cross(xVers);
300
301 G4double xDir = std::sqrt(1. - cosTheta*cosT
302 G4double yDir = xDir;
303 xDir *= std::cos(phi);
304 yDir *= std::sin(phi);
305
306 G4ThreeVector zPrimeVers((xDir*xVers + yDir*
307
308 fParticleChangeForGamma->ProposeMomentumDire
309
310 fParticleChangeForGamma->SetProposedKineticE
311
312 }
313
314 //....oooOO0OOooo........oooOO0OOooo........oo
315
316 G4double G4DNAChampionElasticModel::Theta(G4do
317 G4do
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
333 auto t1 = t2 - 1;
334
335 auto e12 = std::upper_bound(eVecm[(*t1)].beg
336
337
338 auto e11 = e12 - 1;
339
340 auto e22 = std::upper_bound(eVecm[(*t2)].beg
341
342
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][valueE
353 xs12 = eDiffCrossSectionData[valueT1][valueE
354 xs21 = eDiffCrossSectionData[valueT2][valueE
355 xs22 = eDiffCrossSectionData[valueT2][valueE
356
357 if (xs11 == 0 && xs12 == 0 && xs21 == 0 && x
358
359 theta = QuadInterpolator(valueE11, valueE12,
360 xs21, xs22, valueT1
361
362 return theta;
363 }
364
365 //....oooOO0OOooo........oooOO0OOooo........oo
366
367 G4double G4DNAChampionElasticModel::LinLogInte
368
369
370
371
372 {
373 G4double d1 = std::log(xs1);
374 G4double d2 = std::log(xs2);
375 G4double value = G4Exp(d1 + (d2 - d1) * (e -
376 return value;
377 }
378
379 //....oooOO0OOooo........oooOO0OOooo........oo
380
381 G4double G4DNAChampionElasticModel::LinLinInte
382
383
384
385
386 {
387 G4double d1 = xs1;
388 G4double d2 = xs2;
389 G4double value = (d1 + (d2 - d1) * (e - e1)
390 return value;
391 }
392
393 //....oooOO0OOooo........oooOO0OOooo........oo
394
395 G4double G4DNAChampionElasticModel::LogLogInte
396
397
398
399
400 {
401 G4double a = (std::log10(xs2) - std::log10(x
402 / (std::log10(e2) - std::log10(e1));
403 G4double b = std::log10(xs2) - a * std::log1
404 G4double sigma = a * std::log10(e) + b;
405 G4double value = (std::pow(10., sigma));
406 return value;
407 }
408
409 //....oooOO0OOooo........oooOO0OOooo........oo
410
411 G4double G4DNAChampionElasticModel::QuadInterp
412
413
414
415
416
417
418
419
420
421
422
423 {
424 // Log-Log
425 /*
426 G4double interpolatedvalue1 = LogLogInterpo
427 G4double interpolatedvalue2 = LogLogInterpo
428 G4double value = LogLogInterpolate(t1, t2,
429
430
431 // Lin-Log
432 G4double interpolatedvalue1 = LinLogInterpo
433 G4double interpolatedvalue2 = LinLogInterpo
434 G4double value = LinLogInterpolate(t1, t2,
435 */
436
437 // Lin-Lin
438 G4double interpolatedvalue1 = LinLinInterpol
439 G4double interpolatedvalue2 = LinLinInterpol
440 G4double value = LinLinInterpolate(t1, t2, t
441 interpola
442
443 return value;
444 }
445
446 //....oooOO0OOooo........oooOO0OOooo........oo
447
448 G4double G4DNAChampionElasticModel::RandomizeC
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........oo
465
466 void G4DNAChampionElasticModel::SetKillBelowTh
467 {
468 G4ExceptionDescription errMsg;
469 errMsg << "The method G4DNAChampionElasticMo
470
471 G4Exception("G4DNAChampionElasticModel::SetK
472 "deprecated",
473 JustWarning,
474 errMsg);
475 }
476