| 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 // Created on 2016/01/18
27 //
28 // Authors: D. Sakata, W.G. Shin, S. Incerti
29 //
30 // Based on a recent release of the ELSEPA code
31 // developed and provided kindly by F. Salvat et al.
32 // See
33 // Computer Physics Communications, 165(2), 157-190. (2005)
34 // http://dx.doi.org/10.1016/j.cpc.2004.09.006
35 //
36
37 #include "G4DNAELSEPAElasticModel.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4DNAMolecularMaterial.hh"
41 #include "G4EmParameters.hh"
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45 using namespace std;
46
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49 G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel(const G4ParticleDefinition*,
50 const G4String& nam) :
51 G4VEmModel(nam)
52 {
53 verboseLevel = 0;
54
55 G4ProductionCutsTable* theCoupleTable =
56 G4ProductionCutsTable::GetProductionCutsTable();
57 auto numOfCouples = (G4int)theCoupleTable->GetTableSize();
58
59 fpBaseWater = G4Material::GetMaterial("G4_WATER");
60
61 for(G4int i=0; i<numOfCouples; ++i)
62 {
63 const G4MaterialCutsCouple* couple =
64 theCoupleTable->GetMaterialCutsCouple(i);
65
66 const G4Material* material = couple->GetMaterial()->GetBaseMaterial();
67 if(!material) material = couple->GetMaterial();
68
69 auto nelm = (G4int)material->GetNumberOfElements();
70
71 if(nelm==1)
72 {// Protection: only for single element
73 G4int Z = 79;
74 const G4ElementVector* theElementVector = material->GetElementVector();
75 Z = G4lrint((*theElementVector)[0]->GetZ());
76 // Protection: only for GOLD
77 if (Z==79){
78 fkillBelowEnergy_Au = 10. * eV; // Kills e- tracking
79 flowEnergyLimit = 0 * eV; // Must stay at zero for killing
80 fhighEnergyLimit = 1 * GeV; // Default
81 SetLowEnergyLimit (flowEnergyLimit);
82 SetHighEnergyLimit(fhighEnergyLimit);
83 }else{
84 //continue;
85 }
86 }else{// Protection: H2O only is available
87 if(material==fpBaseWater){
88 flowEnergyLimit = 10. * eV;
89 fhighEnergyLimit = 1 * MeV;
90 SetLowEnergyLimit (flowEnergyLimit);
91 SetHighEnergyLimit(fhighEnergyLimit);
92 }else{
93 //continue;
94 }
95 }
96
97 if (verboseLevel > 0)
98 {
99 G4cout << "ELSEPA Elastic model is constructed for "
100 << material->GetName() << G4endl
101 << "Energy range: "
102 << flowEnergyLimit / eV << " eV - "
103 << fhighEnergyLimit / MeV << " MeV"
104 << G4endl;
105 }
106 }
107
108 fParticleChangeForGamma = nullptr;
109 fpMolDensity = nullptr;
110
111 fpData_Au=nullptr;
112 fpData_H2O=nullptr;
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116
117 G4DNAELSEPAElasticModel::~G4DNAELSEPAElasticModel()
118 {
119 delete fpData_Au;
120 delete fpData_H2O;
121
122 eEdummyVec_Au.clear();
123 eEdummyVec_H2O.clear();
124 eCum_Au.clear();
125 eCum_H2O.clear();
126 fAngleData_Au.clear();
127 fAngleData_H2O.clear();
128 }
129
130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132 void G4DNAELSEPAElasticModel::Initialise(const G4ParticleDefinition* particle,
133 const G4DataVector& )
134 {
135 if (verboseLevel > 3)
136 G4cout << "Calling G4DNAELSEPAElasticModel::Initialise()" << G4endl;
137
138 if (isInitialised) {return;}
139
140 if(particle->GetParticleName() != "e-")
141 {
142 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0001",
143 FatalException,"Model not applicable to particle type.");
144 return;
145 }
146
147 G4ProductionCutsTable* theCoupleTable =
148 G4ProductionCutsTable::GetProductionCutsTable();
149 auto numOfCouples = (G4int)theCoupleTable->GetTableSize();
150
151 // UNIT OF TCS
152 G4double scaleFactor = 1.*cm*cm;
153
154 fpData_Au=nullptr;
155 fpData_H2O=nullptr;
156 fpBaseWater = G4Material::GetMaterial("G4_WATER");
157
158 for(G4int i=0; i<numOfCouples; ++i)
159 {
160 const G4MaterialCutsCouple* couple =
161 theCoupleTable->GetMaterialCutsCouple(i);
162 const G4Material* material = couple->GetMaterial()->GetBaseMaterial();
163 if(!material) material = couple->GetMaterial();
164
165 auto nelm = (G4int)material->GetNumberOfElements();
166 if (nelm==1){// Protection: only for single element
167 const G4ElementVector* theElementVector = material->GetElementVector();
168 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
169 if (Z!=79)// Protection: only for GOLD
170 {
171 continue;
172 }
173
174 if (Z>0)
175 {
176 G4String fileZElectron("dna/sigma_elastic_e_elsepa_Z");
177 std::ostringstream oss;
178 oss.str("");
179 oss.clear(stringstream::goodbit);
180 oss << Z;
181 fileZElectron += oss.str()+"_muffintin";
182
183 fpData_Au = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
184 eV,
185 scaleFactor );
186 fpData_Au->LoadData(fileZElectron);
187
188 std::ostringstream eFullFileNameZ;
189 const char *path = G4EmParameters::Instance()->GetDirLEDATA();
190
191 if (path == nullptr)
192 {
193 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0002",
194 FatalException,"G4LEDATA environment variable not set.");
195 return;
196 }
197
198 eFullFileNameZ.str("");
199 eFullFileNameZ.clear(stringstream::goodbit);
200
201 eFullFileNameZ
202 << path
203 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_Z"
204 << Z << "_muffintin.dat";
205
206 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
207
208 if (!eDiffCrossSectionZ)
209 {
210 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0003",
211 FatalException,"Missing data file for cumulated DCS");
212 return;
213 }
214
215 eEdummyVec_Au.clear();
216 eCum_Au.clear();
217 fAngleData_Au.clear();
218
219 eEdummyVec_Au.push_back(0.);
220 do
221 {
222 G4double eDummy;
223 G4double cumDummy;
224 eDiffCrossSectionZ>>eDummy>>cumDummy;
225 if (eDummy != eEdummyVec_Au.back())
226 {
227 eEdummyVec_Au.push_back(eDummy);
228 eCum_Au[eDummy].push_back(0.);
229 }
230 eDiffCrossSectionZ>>fAngleData_Au[eDummy][cumDummy];
231 if (cumDummy != eCum_Au[eDummy].back())
232 {
233 eCum_Au[eDummy].push_back(cumDummy);
234 }
235 }while(!eDiffCrossSectionZ.eof());
236 }
237
238 }else{// Protection: H2O only is available
239 if(material == fpBaseWater && !fpData_H2O){
240 if (LowEnergyLimit() < 10*eV)
241 {
242 G4cout<<"G4DNAELSEPAElasticModel: low energy limit increased from "
243 << LowEnergyLimit()/eV << " eV to " << 10 << " eV"
244 << G4endl;
245 SetLowEnergyLimit(10.*eV);
246 }
247
248 if (HighEnergyLimit() > 1.*MeV)
249 {
250 G4cout<<"G4DNAELSEPAElasticModel: high energy limit decreased from "
251 << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
252 << G4endl;
253 SetHighEnergyLimit(1.*MeV);
254 }
255
256 G4String fileZElectron("dna/sigma_elastic_e_elsepa_muffin");
257
258 fpData_H2O = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
259 eV,
260 scaleFactor );
261 fpData_H2O->LoadData(fileZElectron);
262
263 std::ostringstream eFullFileNameZ;
264
265 const char *path = G4EmParameters::Instance()->GetDirLEDATA();
266 if (path == nullptr)
267 {
268 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0004",
269 FatalException,"G4LEDATA environment variable not set.");
270 return;
271 }
272
273 eFullFileNameZ.str("");
274 eFullFileNameZ.clear(stringstream::goodbit);
275
276 eFullFileNameZ
277 << path
278 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_muffin.dat";
279
280 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
281
282 if (!eDiffCrossSectionZ)
283 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0005",
284 FatalException,
285 "Missing data file for cumulated DCS");
286
287 eEdummyVec_H2O.clear();
288 eCum_H2O.clear();
289 fAngleData_H2O.clear();
290
291 eEdummyVec_H2O.push_back(0.);
292
293 do
294 {
295 G4double eDummy;
296 G4double cumDummy;
297 eDiffCrossSectionZ>>eDummy>>cumDummy;
298 if (eDummy != eEdummyVec_H2O.back())
299 {
300 eEdummyVec_H2O.push_back(eDummy);
301 eCum_H2O[eDummy].push_back(0.);
302 }
303 eDiffCrossSectionZ>>fAngleData_H2O[eDummy][cumDummy];
304 if (cumDummy != eCum_H2O[eDummy].back()){
305 eCum_H2O[eDummy].push_back(cumDummy);
306 }
307 }while(!eDiffCrossSectionZ.eof());
308 }
309 }
310 if (verboseLevel > 2)
311 G4cout << "Loaded cross section files of ELSEPA Elastic model for"
312 << material->GetName() << G4endl;
313
314 if( verboseLevel>0 )
315 {
316 G4cout << "ELSEPA elastic model is initialized " << G4endl
317 << "Energy range: "
318 << LowEnergyLimit() / eV << " eV - "
319 << HighEnergyLimit()/ MeV << " MeV"
320 << G4endl;
321 }
322 } // Loop on couples
323
324
325 fParticleChangeForGamma = GetParticleChangeForGamma();
326
327 fpMolDensity =
328 G4DNAMolecularMaterial::Instance()->
329 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
330
331 isInitialised = true;
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
335
336 G4double G4DNAELSEPAElasticModel::CrossSectionPerVolume
337 (const G4Material* material,
338 const G4ParticleDefinition* particle,
339 G4double ekin,
340 G4double,
341 G4double)
342 {
343
344 if (verboseLevel > 3)
345 {
346 G4cout <<
347 "Calling CrossSectionPerVolume() of G4DNAELSEPAElasticModel"
348 << G4endl;
349 }
350
351 G4double atomicNDensity=0.0;
352 G4double sigma=0;
353
354 std::size_t nelm = material->GetNumberOfElements();
355 if (nelm==1) // Protection: only for single element
356 {
357 // Protection: only for GOLD
358 if (material->GetZ()!=79) return 0.0;
359
360 const G4ElementVector* theElementVector = material->GetElementVector();
361 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
362
363 const G4String& particleName = particle->GetParticleName();
364 atomicNDensity = material->GetAtomicNumDensityVector()[0];
365 if(atomicNDensity!= 0.0)
366 {
367 if (ekin < fhighEnergyLimit)
368 {
369 if (ekin < fkillBelowEnergy_Au) return DBL_MAX;
370
371 if (ekin < 10*eV) sigma = fpData_Au->FindValue(10*eV);
372 else sigma = fpData_Au->FindValue(ekin);
373 }
374 }
375 if (verboseLevel > 2)
376 {
377 G4cout << "__________________________________" << G4endl;
378 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
379 G4cout << "=== Material is made of one element with Z =" << Z << G4endl;
380 G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : "
381 << particleName << G4endl;
382 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^2)"
383 << sigma/cm/cm << G4endl;
384 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^-1)="
385 << sigma*atomicNDensity/(1./cm) << G4endl;
386 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
387 }
388 }
389 else
390 {
391 atomicNDensity = (*fpMolDensity)[material->GetIndex()];
392 if(atomicNDensity!= 0.0)
393 {
394 if (ekin < HighEnergyLimit() && ekin >= LowEnergyLimit())
395 {
396 sigma = fpData_H2O->FindValue(ekin);
397 }
398 }
399 if (verboseLevel > 2)
400 {
401 G4cout << "__________________________________" << G4endl;
402 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
403 G4cout << "=== Kinetic energy(eV)=" << ekin/eV
404 << " particle : " << particle->GetParticleName() << G4endl;
405 G4cout << "=== Cross section per water molecule (cm^2)="
406 << sigma/cm/cm << G4endl;
407 G4cout << "=== Cross section per water molecule (cm^-1)="
408 << sigma*atomicNDensity/(1./cm) << G4endl;
409 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
410 }
411 }
412
413 return sigma*atomicNDensity;
414 }
415
416 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417
418 void G4DNAELSEPAElasticModel::SampleSecondaries(
419 std::vector<G4DynamicParticle*>*,
420 const G4MaterialCutsCouple* couple,
421 const G4DynamicParticle* aDynamicElectron,
422 G4double,
423 G4double)
424 {
425
426 if (verboseLevel > 3){
427 G4cout <<
428 "Calling SampleSecondaries() of G4DNAELSEPAElasticModel"
429 << G4endl;
430 }
431
432 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
433
434 const G4Material* material = couple->GetMaterial()->GetBaseMaterial();
435 if(!material) material = couple->GetMaterial();
436
437 std::size_t nelm = material->GetNumberOfElements();
438 if (nelm==1) // Protection: only for single element
439 {
440 const G4ElementVector* theElementVector = material->GetElementVector();
441 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
442 if (Z!=79) return;
443 if (electronEnergy0 < fkillBelowEnergy_Au)
444 {
445 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
446 fParticleChangeForGamma->ProposeMomentumDirection(G4ThreeVector(0,0,0));
447 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
448 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
449 return;
450 }
451
452 if(electronEnergy0>= fkillBelowEnergy_Au && electronEnergy0 < fhighEnergyLimit)
453 {
454 G4double cosTheta = 0;
455 if (electronEnergy0>=10*eV)
456 {
457 cosTheta = RandomizeCosTheta(Z,electronEnergy0);
458 }
459 else
460 {
461 cosTheta = RandomizeCosTheta(Z,10*eV);
462 }
463
464 G4double phi = 2. * CLHEP::pi * G4UniformRand();
465
466 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
467 G4ThreeVector xVers = zVers.orthogonal();
468 G4ThreeVector yVers = zVers.cross(xVers);
469
470 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
471 G4double yDir = xDir;
472 xDir *= std::cos(phi);
473 yDir *= std::sin(phi);
474
475 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
476 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
477 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
478
479 }
480 }
481 else
482 {
483 if(material == fpBaseWater)
484 {
485 //The data for water is stored as Z=0
486 G4double cosTheta = RandomizeCosTheta(0,electronEnergy0);
487
488 G4double phi = 2. * pi * G4UniformRand();
489
490 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
491 G4ThreeVector xVers = zVers.orthogonal();
492 G4ThreeVector yVers = zVers.cross(xVers);
493
494 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
495 G4double yDir = xDir;
496 xDir *= std::cos(phi);
497 yDir *= std::sin(phi);
498
499 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
500 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
501 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
502 }
503 }
504 }
505
506 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
507
508 G4double G4DNAELSEPAElasticModel::Theta(G4int Z,
509 G4ParticleDefinition * particleDefinition,
510 G4double k,
511 G4double integrDiff)
512 {
513
514 G4double theta = 0.;
515 G4double valueE1 = 0.;
516 G4double valueE2 = 0.;
517 G4double valuecum21 = 0.;
518 G4double valuecum22 = 0.;
519 G4double valuecum12 = 0.;
520 G4double valuecum11 = 0.;
521 G4double a11 = 0.;
522 G4double a12 = 0.;
523 G4double a21 = 0.;
524 G4double a22 = 0.;
525
526 if (particleDefinition == G4Electron::ElectronDefinition())
527 {
528
529 std::vector<G4double>::iterator e2;
530 if(Z==0){
531 e2 = std::upper_bound(eEdummyVec_H2O.begin(),
532 eEdummyVec_H2O.end(), k);
533 }else if (Z==79){
534 e2 = std::upper_bound(eEdummyVec_Au.begin(),
535 eEdummyVec_Au.end(), k);
536 }
537
538 auto e1 = e2 - 1;
539
540 std::vector<G4double>::iterator cum12;
541 if(Z==0){
542 cum12 = std::upper_bound(eCum_H2O[(*e1)].begin(),
543 eCum_H2O[(*e1)].end(),integrDiff);
544 }else if (Z==79){
545 cum12 = std::upper_bound(eCum_Au[(*e1)].begin(),
546 eCum_Au[(*e1)].end(),integrDiff);
547 }
548
549 auto cum11 = cum12 - 1;
550
551 //std::vector<G4double>::iterator cum22
552 // = std::upper_bound(eCumZ[Z][(*e2)].begin(),
553 // eCumZ[Z][(*e2)].end(),integrDiff);
554 std::vector<G4double>::iterator cum22;
555 if(Z==0){
556 cum22 = std::upper_bound(eCum_H2O[(*e2)].begin(),
557 eCum_H2O[(*e2)].end(),integrDiff);
558 }else if(Z==79){
559 cum22 = std::upper_bound(eCum_Au[(*e2)].begin(),
560 eCum_Au[(*e2)].end(),integrDiff);
561 }
562
563 auto cum21 = cum22 - 1;
564
565 valueE1 = *e1;
566 valueE2 = *e2;
567 valuecum11 = *cum11;
568 valuecum12 = *cum12;
569 valuecum21 = *cum21;
570 valuecum22 = *cum22;
571
572 if(Z==0){
573 a11 = fAngleData_H2O[valueE1][valuecum11];
574 a12 = fAngleData_H2O[valueE1][valuecum12];
575 a21 = fAngleData_H2O[valueE2][valuecum21];
576 a22 = fAngleData_H2O[valueE2][valuecum22];
577 }else if (Z==79){
578 a11 = fAngleData_Au[valueE1][valuecum11];
579 a12 = fAngleData_Au[valueE1][valuecum12];
580 a21 = fAngleData_Au[valueE2][valuecum21];
581 a22 = fAngleData_Au[valueE2][valuecum22];
582 }
583
584 }
585
586 if (a11 == 0 && a12 == 0 && a21 == 0 && a22 == 0) return (0.);
587
588 theta = QuadInterpolator(valuecum11, valuecum12, valuecum21, valuecum22,
589 a11, a12,a21, a22, valueE1, valueE2, k, integrDiff);
590 return theta;
591 }
592
593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
594 //
595 G4double G4DNAELSEPAElasticModel::LogLinInterpolate(G4double e1,
596 G4double e2,
597 G4double e,
598 G4double xs1,
599 G4double xs2)
600 {
601 G4double value=0.;
602 if(e1!=0){
603 G4double a = std::log10(e) - std::log10(e1);
604 G4double b = std::log10(e2) - std::log10(e);
605 value = xs1 + a/(a+b)*(xs2-xs1);
606 }
607 else{
608 G4double d1 = xs1;
609 G4double d2 = xs2;
610 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
611 }
612
613 return value;
614 }
615
616 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
617
618 G4double G4DNAELSEPAElasticModel::LinLogInterpolate(G4double e1,
619 G4double e2,
620 G4double e,
621 G4double xs1,
622 G4double xs2)
623 {
624 G4double d1 = std::log10(xs1);
625 G4double d2 = std::log10(xs2);
626 G4double value = std::pow(10,(d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
627 return value;
628 }
629
630 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
631
632 G4double G4DNAELSEPAElasticModel::LinLinInterpolate(G4double e1,
633 G4double e2,
634 G4double e,
635 G4double xs1,
636 G4double xs2)
637 {
638 G4double d1 = xs1;
639 G4double d2 = xs2;
640 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
641 return value;
642 }
643
644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645
646 G4double G4DNAELSEPAElasticModel::LogLogInterpolate(G4double e1,
647 G4double e2,
648 G4double e,
649 G4double xs1,
650 G4double xs2)
651 {
652 G4double a = (std::log10(xs2) - std::log10(xs1))
653 / (std::log10(e2) - std::log10(e1));
654 G4double b = std::log10(xs2) - a * std::log10(e2);
655 G4double sigma = a * std::log10(e) + b;
656 G4double value = (std::pow(10., sigma));
657 return value;
658 }
659
660 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
661
662 G4double G4DNAELSEPAElasticModel::QuadInterpolator(
663 G4double cum11,
664 G4double cum12,
665 G4double cum21,
666 G4double cum22,
667 G4double a11,
668 G4double a12,
669 G4double a21,
670 G4double a22,
671 G4double e1,
672 G4double e2,
673 G4double t,
674 G4double cum)
675 {
676 G4double value=0;
677 G4double interpolatedvalue1=0;
678 G4double interpolatedvalue2=0;
679
680 if(cum11!=0){
681 interpolatedvalue1 = LinLogInterpolate(cum11, cum12, cum, a11, a12);
682 }
683 else{
684 interpolatedvalue1 = LinLinInterpolate(cum11, cum12, cum, a11, a12);
685 }
686 if(cum21!=0){
687 interpolatedvalue2 = LinLogInterpolate(cum21, cum22, cum, a21, a22);
688 }
689 else{
690 interpolatedvalue2 = LinLinInterpolate(cum21, cum22, cum, a21, a22);
691 }
692
693 value = LogLinInterpolate(e1,e2,t,interpolatedvalue1,interpolatedvalue2);
694
695 return value;
696 }
697
698 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
699
700 G4double G4DNAELSEPAElasticModel::RandomizeCosTheta(G4int Z, G4double k)
701 {
702
703 G4double integrdiff = 0.;
704 G4double uniformRand = G4UniformRand();
705 integrdiff = uniformRand;
706
707 G4double theta = 0.;
708 G4double cosTheta = 0.;
709 theta = Theta(Z, G4Electron::ElectronDefinition(), k / eV, integrdiff);
710
711 cosTheta = std::cos(theta * CLHEP::pi / 180.);
712
713 return cosTheta;
714 }
715
716 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
717
718 void G4DNAELSEPAElasticModel::SetKillBelowThreshold(G4double threshold)
719 {
720 fkillBelowEnergy_Au = threshold;
721
722 if (threshold < 10 * eV)
723 {
724 G4cout<< "*** WARNING : the G4DNAELSEPAElasticModel model is not "
725 "defined below 10 eV !" << G4endl;
726 }
727 }
728