Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAELSEPAElasticModel.cc
Go to the documentation of this file.
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
45using namespace std;
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel(const G4ParticleDefinition*,
50const G4String& nam) :
51G4VEmModel(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
117G4DNAELSEPAElasticModel::~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
132void G4DNAELSEPAElasticModel::Initialise(const G4ParticleDefinition* particle,
133const 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
336G4double 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
418void 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
508G4double 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//
595G4double G4DNAELSEPAElasticModel::LogLinInterpolate(G4double e1,
596G4double e2,
597G4double e,
598G4double xs1,
599G4double 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
618G4double G4DNAELSEPAElasticModel::LinLogInterpolate(G4double e1,
619G4double e2,
620G4double e,
621G4double xs1,
622G4double 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
632G4double G4DNAELSEPAElasticModel::LinLinInterpolate(G4double e1,
633G4double e2,
634G4double e,
635G4double xs1,
636G4double 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
646G4double G4DNAELSEPAElasticModel::LogLogInterpolate(G4double e1,
647G4double e2,
648G4double e,
649G4double xs1,
650G4double 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
662G4double G4DNAELSEPAElasticModel::QuadInterpolator(
663G4double cum11,
664G4double cum12,
665G4double cum21,
666G4double cum22,
667G4double a11,
668G4double a12,
669G4double a21,
670G4double a22,
671G4double e1,
672G4double e2,
673G4double t,
674G4double 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
700G4double 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
718void 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}
G4ElementVector
std::vector< const G4Element * > G4ElementVector
Definition G4ElementVector.hh:33
FatalException
@ FatalException
Definition G4ExceptionSeverity.hh:67
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition G4Exception.cc:59
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition G4ThreeVector.hh:36
fStopAndKill
@ fStopAndKill
Definition G4TrackStatus.hh:46
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52
CLHEP::Hep3Vector::unit
Hep3Vector unit() const
CLHEP::Hep3Vector::orthogonal
Hep3Vector orthogonal() const
CLHEP::Hep3Vector::cross
Hep3Vector cross(const Hep3Vector &) const
G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel
G4DNAELSEPAElasticModel(const G4ParticleDefinition *particle=nullptr, const G4String &nam="DNAELSEPAElasticModel")
Definition G4DNAELSEPAElasticModel.cc:49
G4DNAELSEPAElasticModel::CrossSectionPerVolume
G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *particle, G4double ekin, G4double emin, G4double emax) override
Definition G4DNAELSEPAElasticModel.cc:337
G4DNAELSEPAElasticModel::Initialise
void Initialise(const G4ParticleDefinition *particle, const G4DataVector &) override
Definition G4DNAELSEPAElasticModel.cc:132
G4DNAELSEPAElasticModel::SampleSecondaries
void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
Definition G4DNAELSEPAElasticModel.cc:418
G4DNAELSEPAElasticModel::SetKillBelowThreshold
void SetKillBelowThreshold(G4double threshold)
Definition G4DNAELSEPAElasticModel.cc:718
G4DNAELSEPAElasticModel::fParticleChangeForGamma
G4ParticleChangeForGamma * fParticleChangeForGamma
Definition G4DNAELSEPAElasticModel.hh:86
G4DNAMolecularMaterial::Instance
static G4DNAMolecularMaterial * Instance()
Definition G4DNAMolecularMaterial.cc:80
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4Electron::ElectronDefinition
static G4Electron * ElectronDefinition()
Definition G4Electron.cc:86
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition G4EmParameters.cc:70
G4EmParameters::GetDirLEDATA
const G4String & GetDirLEDATA() const
Definition G4EmParameters.cc:1363
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition G4MaterialCutsCouple.hh:166
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition G4Material.hh:183
G4Material::GetBaseMaterial
const G4Material * GetBaseMaterial() const
Definition G4Material.hh:218
G4Material::GetZ
G4double GetZ() const
Definition G4Material.cc:710
G4Material::GetAtomicNumDensityVector
const G4double * GetAtomicNumDensityVector() const
Definition G4Material.hh:202
G4Material::GetIndex
std::size_t GetIndex() const
Definition G4Material.hh:239
G4Material::GetNumberOfElements
std::size_t GetNumberOfElements() const
Definition G4Material.hh:180
G4Material::GetName
const G4String & GetName() const
Definition G4Material.hh:171
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition G4Material.cc:663
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition G4ParticleDefinition.hh:96
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition G4ProductionCutsTable.hh:256
G4ProductionCutsTable::GetTableSize
std::size_t GetTableSize() const
Definition G4ProductionCutsTable.hh:249
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition G4ProductionCutsTable.cc:58
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition G4VEmModel.hh:731
G4VEmModel::GetParticleChangeForGamma
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition G4VEmModel.cc:121
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition G4VEmModel.hh:640
G4VEmModel::HighEnergyLimit
G4double HighEnergyLimit() const
Definition G4VEmModel.hh:633
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition G4VEmModel.hh:738
G4VEmModel::G4VEmModel
G4VEmModel(const G4String &nam)
Definition G4VEmModel.cc:67
G4lrint
int G4lrint(double ad)
Definition templates.hh:134
DBL_MAX
#define DBL_MAX
Definition templates.hh:62