| 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 // 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 cod
31 // developed and provided kindly by F. Salvat
32 // See
33 // Computer Physics Communications, 165(2), 15
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........oo
44
45 using namespace std;
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4DNAELSEPAElasticModel::G4DNAELSEPAElasticMod
50 const G4String& nam) :
51 G4VEmModel(nam)
52 {
53 verboseLevel = 0;
54
55 G4ProductionCutsTable* theCoupleTable =
56 G4ProductionCutsTable::GetProductionCutsTabl
57 auto numOfCouples = (G4int)theCoupleTable->
58
59 fpBaseWater = G4Material::GetMaterial("G4_WA
60
61 for(G4int i=0; i<numOfCouples; ++i)
62 {
63 const G4MaterialCutsCouple* couple =
64 theCoupleTable->GetMaterialCutsCouple
65
66 const G4Material* material = couple->GetMa
67 if(!material) material = couple->GetMateri
68
69 auto nelm = (G4int)material->GetNumberOfE
70
71 if(nelm==1)
72 {// Protection: only for single element
73 G4int Z = 79;
74 const G4ElementVector* theElementVector
75 Z = G4lrint((*theElementVector)[0]->Get
76 // Protection: only for GOLD
77 if (Z==79){
78 fkillBelowEnergy_Au = 10. * eV; // Ki
79 flowEnergyLimit = 0 * eV; // Must
80 fhighEnergyLimit = 1 * GeV; // Defau
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 const
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........oo
116
117 G4DNAELSEPAElasticModel::~G4DNAELSEPAElasticMo
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........oo
131
132 void G4DNAELSEPAElasticModel::Initialise(const
133 const G4DataVector& )
134 {
135 if (verboseLevel > 3)
136 G4cout << "Calling G4DNAELSEPAElasticModel::
137
138 if (isInitialised) {return;}
139
140 if(particle->GetParticleName() != "e-")
141 {
142 G4Exception("G4DNAELSEPAElasticModel::Init
143 FatalException,"Model not applicable to
144 return;
145 }
146
147 G4ProductionCutsTable* theCoupleTable =
148 G4ProductionCutsTable::GetProductionCutsTabl
149 auto numOfCouples = (G4int)theCoupleTable->
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_WA
157
158 for(G4int i=0; i<numOfCouples; ++i)
159 {
160 const G4MaterialCutsCouple* couple =
161 theCoupleTable->GetMaterialCutsCouple
162 const G4Material* material = couple->GetMa
163 if(!material) material = couple->GetMateri
164
165 auto nelm = (G4int)material->GetNumberOfE
166 if (nelm==1){// Protection: only for singl
167 const G4ElementVector* theElementVector
168 G4int Z = G4lrint((*theElementVector)[0
169 if (Z!=79)// Protection: only for GOLD
170 {
171 continue;
172 }
173
174 if (Z>0)
175 {
176 G4String fileZElectron("dna/sigma_elas
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 G4DNACrossSectionDataS
184
185
186 fpData_Au->LoadData(fileZElectron);
187
188 std::ostringstream eFullFileNameZ;
189 const char *path = G4EmParameters::Ins
190
191 if (path == nullptr)
192 {
193 G4Exception("G4DNAELSEPAElasticModel
194 FatalException,"G4LEDATA environme
195 return;
196 }
197
198 eFullFileNameZ.str("");
199 eFullFileNameZ.clear(stringstream::goo
200
201 eFullFileNameZ
202 << path
203 << "/dna/sigmadiff_cumulated_elastic
204 << Z << "_muffintin.dat";
205
206 std::ifstream eDiffCrossSectionZ(eFull
207
208 if (!eDiffCrossSectionZ)
209 {
210 G4Exception("G4DNAELSEPAElasticModel
211 FatalException,"Missing data file
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[eD
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_H2
240 if (LowEnergyLimit() < 10*eV)
241 {
242 G4cout<<"G4DNAELSEPAElasticModel: lo
243 << LowEnergyLimit()/eV << " eV
244 << G4endl;
245 SetLowEnergyLimit(10.*eV);
246 }
247
248 if (HighEnergyLimit() > 1.*MeV)
249 {
250 G4cout<<"G4DNAELSEPAElasticModel: hi
251 << HighEnergyLimit()/MeV << "
252 << G4endl;
253 SetHighEnergyLimit(1.*MeV);
254 }
255
256 G4String fileZElectron("dna/sigma_elas
257
258 fpData_H2O = new G4DNACrossSectionData
259
260
261 fpData_H2O->LoadData(fileZElectron);
262
263 std::ostringstream eFullFileNameZ;
264
265 const char *path = G4EmParameters::Ins
266 if (path == nullptr)
267 {
268 G4Exception("G4DNAELSEPAElasticModel
269 FatalException,"G4LEDATA environme
270 return;
271 }
272
273 eFullFileNameZ.str("");
274 eFullFileNameZ.clear(stringstream::goo
275
276 eFullFileNameZ
277 << path
278 << "/dna/sigmadiff_cumulated_elasti
279
280 std::ifstream eDiffCrossSectionZ(eFull
281
282 if (!eDiffCrossSectionZ)
283 G4Exception("G4DNAELSEPAElasticModel:
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[e
304 if (cumDummy != eCum_H2O[eDummy].bac
305 eCum_H2O[eDummy].push_back(cumDumm
306 }
307 }while(!eDiffCrossSectionZ.eof());
308 }
309 }
310 if (verboseLevel > 2)
311 G4cout << "Loaded cross section files of E
312 << material->GetName() << G4endl;
313
314 if( verboseLevel>0 )
315 {
316 G4cout << "ELSEPA elastic model is initi
317 << "Energy range: "
318 << LowEnergyLimit() / eV << " eV - "
319 << HighEnergyLimit()/ MeV << " MeV"
320 << G4endl;
321 }
322 } // Loop on couples
323
324
325 fParticleChangeForGamma = GetParticleChangeF
326
327 fpMolDensity =
328 G4DNAMolecularMaterial::Instance()->
329 GetNumMolPerVolTableFor(G4Material::GetMater
330
331 isInitialised = true;
332 }
333
334 //....oooOO0OOooo........oooOO0OOooo........oo
335
336 G4double G4DNAELSEPAElasticModel::CrossSection
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 G4DNAE
348 << G4endl;
349 }
350
351 G4double atomicNDensity=0.0;
352 G4double sigma=0;
353
354 std::size_t nelm = material->GetNumberOfElem
355 if (nelm==1) // Protection: only for single
356 {
357 // Protection: only for GOLD
358 if (material->GetZ()!=79) return 0.0;
359
360 const G4ElementVector* theElementVector =
361 G4int Z = G4lrint((*theElementVector)[0]->
362
363 const G4String& particleName = particle->G
364 atomicNDensity = material->GetAtomicNumDen
365 if(atomicNDensity!= 0.0)
366 {
367 if (ekin < fhighEnergyLimit)
368 {
369 if (ekin < fkillBelowEnergy_Au) return
370
371 if (ekin < 10*eV) sigma = fpData_Au->F
372 else sigma = fpData_Au->F
373 }
374 }
375 if (verboseLevel > 2)
376 {
377 G4cout << "_____________________________
378 G4cout << "=== G4DNAELSEPAElasticModel -
379 G4cout << "=== Material is made of one e
380 G4cout << "=== Kinetic energy(eV)=" << e
381 << particleName << G4endl;
382 G4cout << "=== Cross section per atom fo
383 << sigma/cm/cm << G4endl;
384 G4cout << "=== Cross section per atom fo
385 << sigma*atomicNDensity/(1./cm) <
386 G4cout << "=== G4DNAELSEPAElasticModel -
387 }
388 }
389 else
390 {
391 atomicNDensity = (*fpMolDensity)[material-
392 if(atomicNDensity!= 0.0)
393 {
394 if (ekin < HighEnergyLimit() && ekin >=
395 {
396 sigma = fpData_H2O->FindValue(ekin);
397 }
398 }
399 if (verboseLevel > 2)
400 {
401 G4cout << "_____________________________
402 G4cout << "=== G4DNAELSEPAElasticModel -
403 G4cout << "=== Kinetic energy(eV)=" << e
404 << " particle : " << particle->Ge
405 G4cout << "=== Cross section per water m
406 << sigma/cm/cm << G4endl;
407 G4cout << "=== Cross section per water m
408 << sigma*atomicNDensity/(1./cm) <
409 G4cout << "=== G4DNAELSEPAElasticModel -
410 }
411 }
412
413 return sigma*atomicNDensity;
414 }
415
416 //....oooOO0OOooo........oooOO0OOooo........oo
417
418 void G4DNAELSEPAElasticModel::SampleSecondarie
419 std::vector<G4DynamicParticle*>*,
420 const G4MaterialCutsCouple* couple,
421 const G4DynamicParticle* aDynamicElectro
422 G4double,
423 G4double)
424 {
425
426 if (verboseLevel > 3){
427 G4cout <<
428 "Calling SampleSecondaries() of G4DNAELSEP
429 << G4endl;
430 }
431
432 G4double electronEnergy0 = aDynamicElectron-
433
434 const G4Material* material = couple->GetMate
435 if(!material) material = couple->GetMaterial
436
437 std::size_t nelm = material->GetNumberOfElem
438 if (nelm==1) // Protection: only for single
439 {
440 const G4ElementVector* theElementVector =
441 G4int Z = G4lrint((*theElementVector)[0]-
442 if (Z!=79) return;
443 if (electronEnergy0 < fkillBelowEnergy_Au)
444 {
445 fParticleChangeForGamma->SetProposedKine
446 fParticleChangeForGamma->ProposeMomentum
447 fParticleChangeForGamma->ProposeTrackSta
448 fParticleChangeForGamma->ProposeLocalEne
449 return;
450 }
451
452 if(electronEnergy0>= fkillBelowEnergy_Au &
453 {
454 G4double cosTheta = 0;
455 if (electronEnergy0>=10*eV)
456 {
457 cosTheta = RandomizeCosTheta(Z,electro
458 }
459 else
460 {
461 cosTheta = RandomizeCosTheta(Z,10*eV);
462 }
463
464 G4double phi = 2. * CLHEP::pi * G4Unifor
465
466 G4ThreeVector zVers = aDynamicElectron->
467 G4ThreeVector xVers = zVers.orthogonal()
468 G4ThreeVector yVers = zVers.cross(xVers)
469
470 G4double xDir = std::sqrt(1. - cosTheta*
471 G4double yDir = xDir;
472 xDir *= std::cos(phi);
473 yDir *= std::sin(phi);
474
475 G4ThreeVector zPrimeVers((xDir*xVers + y
476 fParticleChangeForGamma->ProposeMomentum
477 fParticleChangeForGamma->SetProposedKine
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,
487
488 G4double phi = 2. * pi * G4UniformRand()
489
490 G4ThreeVector zVers = aDynamicElectron->
491 G4ThreeVector xVers = zVers.orthogonal()
492 G4ThreeVector yVers = zVers.cross(xVers)
493
494 G4double xDir = std::sqrt(1. - cosTheta*
495 G4double yDir = xDir;
496 xDir *= std::cos(phi);
497 yDir *= std::sin(phi);
498
499 G4ThreeVector zPrimeVers((xDir*xVers + y
500 fParticleChangeForGamma->ProposeMomentum
501 fParticleChangeForGamma->SetProposedKine
502 }
503 }
504 }
505
506 //....oooOO0OOooo........oooOO0OOooo........oo
507
508 G4double G4DNAELSEPAElasticModel::Theta(G4int
509 G4ParticleDefinition * particleDefini
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::Electro
527 {
528
529 std::vector<G4double>::iterator e2;
530 if(Z==0){
531 e2 = std::upper_bound(eEdummyVec_H2O.begin
532 eEdummyVec_H2O.end()
533 }else if (Z==79){
534 e2 = std::upper_bound(eEdummyVec_Au.begin(
535 eEdummyVec_Au.end(),
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)]
543 eCum_H2O[(*e1)]
544 }else if (Z==79){
545 cum12 = std::upper_bound(eCum_Au[(*e1)].
546 eCum_Au[(*e1)].
547 }
548
549 auto cum11 = cum12 - 1;
550
551 //std::vector<G4double>::iterator cum22
552 // = std::upper_bound(eCumZ[Z][(*e
553 // eCumZ[Z][(*e2)].end(),integrDif
554 std::vector<G4double>::iterator cum22;
555 if(Z==0){
556 cum22 = std::upper_bound(eCum_H2O[(*e2)].
557 eCum_H2O[(*e2)].
558 }else if(Z==79){
559 cum22 = std::upper_bound(eCum_Au[(*e2)].b
560 eCum_Au[(*e2)].e
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 =
587
588 theta = QuadInterpolator(valuecum11, valuecum
589 a11, a12,a21, a22, valueE1, valueE2,
590 return theta;
591 }
592
593 //....oooOO0OOooo........oooOO0OOooo........oo
594 //
595 G4double G4DNAELSEPAElasticModel::LogLinInterp
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 -
611 }
612
613 return value;
614 }
615
616 //....oooOO0OOooo........oooOO0OOooo........oo
617
618 G4double G4DNAELSEPAElasticModel::LinLogInterp
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)
627 return value;
628 }
629
630 //....oooOO0OOooo........oooOO0OOooo........oo
631
632 G4double G4DNAELSEPAElasticModel::LinLinInterp
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) /
641 return value;
642 }
643
644 //....oooOO0OOooo........oooOO0OOooo........oo
645
646 G4double G4DNAELSEPAElasticModel::LogLogInterp
647 G4double e2,
648 G4double e,
649 G4double xs1,
650 G4double xs2)
651 {
652 G4double a = (std::log10(xs2) - std::log10(xs
653 / (std::log10(e2) - std::log10(e1
654 G4double b = std::log10(xs2) - a * std::log10
655 G4double sigma = a * std::log10(e) + b;
656 G4double value = (std::pow(10., sigma));
657 return value;
658 }
659
660 //....oooOO0OOooo........oooOO0OOooo........oo
661
662 G4double G4DNAELSEPAElasticModel::QuadInterpol
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(cu
682 }
683 else{
684 interpolatedvalue1 = LinLinInterpolate(cu
685 }
686 if(cum21!=0){
687 interpolatedvalue2 = LinLogInterpolate(cu
688 }
689 else{
690 interpolatedvalue2 = LinLinInterpolate(cu
691 }
692
693 value = LogLinInterpolate(e1,e2,t,interpola
694
695 return value;
696 }
697
698 //....oooOO0OOooo........oooOO0OOooo........oo
699
700 G4double G4DNAELSEPAElasticModel::RandomizeCos
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::ElectronDefinit
710
711 cosTheta = std::cos(theta * CLHEP::pi / 180.
712
713 return cosTheta;
714 }
715
716 //....oooOO0OOooo........oooOO0OOooo........oo
717
718 void G4DNAELSEPAElasticModel::SetKillBelowThre
719 {
720 fkillBelowEnergy_Au = threshold;
721
722 if (threshold < 10 * eV)
723 {
724 G4cout<< "*** WARNING : the G4DNAELSEPAEla
725 "defined below 10 eV !" << G4endl;
726 }
727 }
728