| 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 "G4DNABornIonisationModel1.hh"
29 #include "G4PhysicalConstants.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4UAtomicDeexcitation.hh"
32 #include "G4LossTableManager.hh"
33 #include "G4DNAChemistryManager.hh"
34 #include "G4DNAMolecularMaterial.hh"
35 #include "G4DNABornAngle.hh"
36 #include "G4DeltaAngle.hh"
37 #include "G4Exp.hh"
38
39 //....oooOO0OOooo........oooOO0OOooo........oo
40
41 using namespace std;
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 G4DNABornIonisationModel1::G4DNABornIonisation
46
47 G4VEmModel(nam)
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 if (verboseLevel > 0)
58 {
59 G4cout << "Born ionisation model is constr
60 }
61
62 // Mark this model as "applicable" for atomi
63 SetDeexcitationFlag(true);
64 fAtomDeexcitation = nullptr;
65 fParticleChangeForGamma = nullptr;
66 fpMolWaterDensity = nullptr;
67
68 // Define default angular generator
69 SetAngularDistribution(new G4DNABornAngle())
70
71 // Selection of computation method
72
73 fasterCode = false;
74
75 // Selection of stationary mode
76
77 statCode = false;
78
79 // Selection of SP scaling
80
81 spScaling = true;
82 }
83
84 //....oooOO0OOooo........oooOO0OOooo........oo
85
86 G4DNABornIonisationModel1::~G4DNABornIonisatio
87 {
88 // Cross section
89
90 std::map<G4String, G4DNACrossSectionDataSet*
91 for (pos = tableData.begin(); pos != tableDa
92 {
93 G4DNACrossSectionDataSet* table = pos->sec
94 delete table;
95 }
96
97 // Final state
98
99 eVecm.clear();
100 pVecm.clear();
101 }
102
103 //....oooOO0OOooo........oooOO0OOooo........oo
104
105 void G4DNABornIonisationModel1::Initialise(con
106 con
107 {
108
109 if (verboseLevel > 3)
110 {
111 G4cout << "Calling G4DNABornIonisationMode
112 }
113
114 // Energy limits
115
116 G4String fileElectron("dna/sigma_ionisation_
117 G4String fileProton("dna/sigma_ionisation_p_
118
119 G4ParticleDefinition* electronDef = G4Electr
120 G4ParticleDefinition* protonDef = G4Proton::
121
122 G4String electron;
123 G4String proton;
124
125 G4double scaleFactor = (1.e-22 / 3.343) * m*
126
127 const char *path = G4FindDataDir("G4LEDATA")
128
129 // *** ELECTRON
130
131 electron = electronDef->GetParticleName();
132
133 tableFile[electron] = fileElectron;
134
135 lowEnergyLimit[electron] = 11. * eV;
136 highEnergyLimit[electron] = 1. * MeV;
137
138 // Cross section
139
140 auto tableE = new G4DNACrossSectionDataSet(
141 tableE->LoadData(fileElectron);
142
143 tableData[electron] = tableE;
144
145 // Final state
146
147 std::ostringstream eFullFileName;
148
149 if (fasterCode) eFullFileName << path << "/d
150 if (!fasterCode) eFullFileName << path << "/
151
152 std::ifstream eDiffCrossSection(eFullFileNam
153
154 if (!eDiffCrossSection)
155 {
156 if (fasterCode) G4Exception("G4DNABornIoni
157 FatalException,"Missing data file:/dna
158
159 if (!fasterCode) G4Exception("G4DNABornIon
160 FatalException,"Missing data file:/dna
161 }
162
163 // Clear the arrays for re-initialization ca
164 // March 25th, 2014 - Vaclav Stepan, Sebasti
165
166 eTdummyVec.clear();
167 pTdummyVec.clear();
168
169 eVecm.clear();
170 pVecm.clear();
171
172 for (G4int j=0; j<5; j++)
173 {
174 eProbaShellMap[j].clear();
175 pProbaShellMap[j].clear();
176
177 eDiffCrossSectionData[j].clear();
178 pDiffCrossSectionData[j].clear();
179
180 eNrjTransfData[j].clear();
181 pNrjTransfData[j].clear();
182 }
183
184 //
185
186 eTdummyVec.push_back(0.);
187 while(!eDiffCrossSection.eof())
188 {
189 G4double tDummy;
190 G4double eDummy;
191 eDiffCrossSection>>tDummy>>eDummy;
192 if (tDummy != eTdummyVec.back()) eTdummyVe
193
194 G4double tmp;
195 for (G4int j=0; j<5; j++)
196 {
197 eDiffCrossSection>> tmp;
198
199 eDiffCrossSectionData[j][tDummy][eDummy]
200
201 if (fasterCode)
202 {
203 eNrjTransfData[j][tDummy][eDiffCrossSe
204 eProbaShellMap[j][tDummy].push_back(eD
205 }
206
207 // SI - only if eof is not reached
208 if (!eDiffCrossSection.eof() && !fasterC
209
210 if (!fasterCode) eVecm[tDummy].push_back
211
212 }
213 }
214
215 // *** PROTON
216
217 proton = protonDef->GetParticleName();
218
219 tableFile[proton] = fileProton;
220
221 lowEnergyLimit[proton] = 500. * keV;
222 highEnergyLimit[proton] = 100. * MeV;
223
224 // Cross section
225
226 auto tableP = new G4DNACrossSectionDataSet(
227 tableP->LoadData(fileProton);
228
229 tableData[proton] = tableP;
230
231 // Final state
232
233 std::ostringstream pFullFileName;
234
235 if (fasterCode) pFullFileName << path << "/d
236
237 if (!fasterCode) pFullFileName << path << "/
238
239 std::ifstream pDiffCrossSection(pFullFileNam
240
241 if (!pDiffCrossSection)
242 {
243 if (fasterCode) G4Exception("G4DNABornIoni
244 FatalException,"Missing data file:/dna
245
246 if (!fasterCode) G4Exception("G4DNABornIon
247 FatalException,"Missing data file:/dna
248 }
249
250 pTdummyVec.push_back(0.);
251 while(!pDiffCrossSection.eof())
252 {
253 G4double tDummy;
254 G4double eDummy;
255 pDiffCrossSection>>tDummy>>eDummy;
256 if (tDummy != pTdummyVec.back()) pTdummyVe
257 for (G4int j=0; j<5; j++)
258 {
259 pDiffCrossSection>>pDiffCrossSectionData
260
261 if (fasterCode)
262 {
263 pNrjTransfData[j][tDummy][pDiffCrossSe
264 pProbaShellMap[j][tDummy].push_back(pD
265 }
266
267 // SI - only if eof is not reached !
268 if (!pDiffCrossSection.eof() && !fasterC
269
270 if (!fasterCode) pVecm[tDummy].push_back
271 }
272 }
273
274 //
275
276 if (particle==electronDef)
277 {
278 SetLowEnergyLimit(lowEnergyLimit[electron]
279 SetHighEnergyLimit(highEnergyLimit[electro
280 }
281
282 if (particle==protonDef)
283 {
284 SetLowEnergyLimit(lowEnergyLimit[proton]);
285 SetHighEnergyLimit(highEnergyLimit[proton]
286 }
287
288 if( verboseLevel>0 )
289 {
290 G4cout << "Born ionisation model is initia
291 << "Energy range: "
292 << LowEnergyLimit() / eV << " eV - "
293 << HighEnergyLimit() / keV << " keV for "
294 << particle->GetParticleName()
295 << G4endl;
296 }
297
298 // Initialize water density pointer
299
300 fpMolWaterDensity = G4DNAMolecularMaterial::
301 GetNumMolPerVolTableFor(G4Material::GetMater
302
303 // AD
304
305 fAtomDeexcitation = G4LossTableManager::Inst
306
307 //
308
309 if (isInitialised)
310 { return;}
311 fParticleChangeForGamma = GetParticleChangeF
312 isInitialised = true;
313 }
314
315 //....oooOO0OOooo........oooOO0OOooo........oo
316
317 G4double G4DNABornIonisationModel1::CrossSecti
318
319
320
321
322 {
323 if (verboseLevel > 3)
324 {
325 G4cout << "Calling CrossSectionPerVolume()
326 << G4endl;
327 }
328
329 if (
330 particleDefinition != G4Proton::ProtonDe
331 &&
332 particleDefinition != G4Electron::Electr
333 )
334
335 return 0;
336
337 // Calculate total cross section for model
338
339 G4double lowLim = 0;
340 G4double highLim = 0;
341 G4double sigma=0;
342
343 G4double waterDensity = (*fpMolWaterDensity)
344
345 const G4String& particleName = particleDefin
346
347 std::map< G4String,G4double,std::less<G4Stri
348 pos1 = lowEnergyLimit.find(particleName);
349 if (pos1 != lowEnergyLimit.end())
350 {
351 lowLim = pos1->second;
352 }
353
354 std::map< G4String,G4double,std::less<G4Stri
355 pos2 = highEnergyLimit.find(particleName);
356 if (pos2 != highEnergyLimit.end())
357 {
358 highLim = pos2->second;
359 }
360
361 if (ekin >= lowLim && ekin <= highLim)
362 {
363 std::map< G4String,G4DNACrossSectionDataSe
364 pos = tableData.find(particleName);
365
366 if (pos != tableData.end())
367 {
368 G4DNACrossSectionDataSet* table = pos->s
369 if (table != nullptr)
370 {
371 sigma = table->FindValue(ekin);
372
373 // ICRU49 electronic SP scaling - ZF,
374
375 if (particleDefinition == G4Proton::Pr
376 {
377 G4double A = 1.39241700556072800000E-
378 G4double B = -8.52610412942622630000E
379 sigma = sigma * G4Exp(A*(ekin/eV)+B);
380 }
381 //
382
383 }
384 }
385 else
386 {
387 G4Exception("G4DNABornIonisationModel1::
388 FatalException,"Model not applicable
389 }
390 }
391
392 if (verboseLevel > 2)
393 {
394 G4cout << "_______________________________
395 G4cout << "G4DNABornIonisationModel1 - XS
396 G4cout << "Kinetic energy(eV)=" << ekin/eV
397 G4cout << "Cross section per water molecul
398 G4cout << "Cross section per water molecul
399 G4cout << "G4DNABornIonisationModel1 - XS
400 }
401
402 return sigma*waterDensity;
403 }
404
405 //....oooOO0OOooo........oooOO0OOooo........oo
406
407 void G4DNABornIonisationModel1::SampleSecondar
408
409
410
411
412 {
413
414 if (verboseLevel > 3)
415 {
416 G4cout << "Calling SampleSecondaries() of
417 << G4endl;
418 }
419
420 G4double lowLim = 0;
421 G4double highLim = 0;
422
423 G4double k = particle->GetKineticEnergy();
424
425 const G4String& particleName = particle->Get
426
427 std::map< G4String,G4double,std::less<G4Stri
428 pos1 = lowEnergyLimit.find(particleName);
429
430 if (pos1 != lowEnergyLimit.end())
431 {
432 lowLim = pos1->second;
433 }
434
435 std::map< G4String,G4double,std::less<G4Stri
436 pos2 = highEnergyLimit.find(particleName);
437
438 if (pos2 != highEnergyLimit.end())
439 {
440 highLim = pos2->second;
441 }
442
443 if (k >= lowLim && k <= highLim)
444 {
445 G4ParticleMomentum primaryDirection = part
446 G4double particleMass = particle->GetDefin
447 G4double totalEnergy = k + particleMass;
448 G4double pSquare = k * (totalEnergy + part
449 G4double totalMomentum = std::sqrt(pSquare
450
451 G4int ionizationShell = 0;
452
453 if (!fasterCode) ionizationShell = RandomS
454
455 // SI: The following protection is necessa
456 // sigmadiff_ionisation_e_born.dat has no
457 // sigmadiff_cumulated_ionisation_e_born.
458 // this is due to the fact that the max a
459 // strictly above this value have non zer
460
461 if (fasterCode)
462 do
463 {
464 ionizationShell = RandomSelect(k,particl
465 } while (k<19*eV && ionizationShell==2 &&
466
467 G4double bindingEnergy = 0;
468 bindingEnergy = waterStructure.IonisationE
469
470 // SI: additional protection if tcs interp
471 if (k<bindingEnergy) return;
472 //
473
474 G4double secondaryKinetic=-1000*eV;
475
476 if (!fasterCode)
477 {
478 secondaryKinetic = RandomizeEjectedElect
479 }
480 else
481 {
482 secondaryKinetic = RandomizeEjectedElect
483 }
484 //
485
486 G4int Z = 8;
487
488 G4ThreeVector deltaDirection =
489 GetAngularDistribution()->SampleDirectionF
490 Z, ionizationShell,
491 couple->GetMaterial());
492
493 if (secondaryKinetic>0)
494 {
495 auto dp = new G4DynamicParticle (G4Elec
496 fvect->push_back(dp);
497 }
498
499 if (particle->GetDefinition() == G4Electro
500 {
501 G4double deltaTotalMomentum = std::sqrt(
502
503 G4double finalPx = totalMomentum*primary
504 G4double finalPy = totalMomentum*primary
505 G4double finalPz = totalMomentum*primary
506 G4double finalMomentum = std::sqrt(final
507 finalPx /= finalMomentum;
508 finalPy /= finalMomentum;
509 finalPz /= finalMomentum;
510
511 G4ThreeVector direction;
512 direction.set(finalPx,finalPy,finalPz);
513
514 fParticleChangeForGamma->ProposeMomentum
515 }
516
517 else fParticleChangeForGamma->ProposeMomen
518
519 // AM: sample deexcitation
520 // here we assume that H_{2}O electronic l
521 // this can be considered true with a roug
522
523 std::size_t secNumberInit = 0;// need to k
524 std::size_t secNumberFinal = 0;// So I'll
525
526 G4double scatteredEnergy = k-bindingEnergy
527
528 // SI: only atomic deexcitation from K she
529 if((fAtomDeexcitation != nullptr) && ioniz
530 {
531 const G4AtomicShell* shell =
532 fAtomDeexcitation->GetAtomicShell(Z, G
533 secNumberInit = fvect->size();
534 fAtomDeexcitation->GenerateParticles(fve
535 secNumberFinal = fvect->size();
536
537 //TEST
538 //G4cout << "ionizationShell=" << ioniza
539 //G4cout << "bindingEnergy=" << bindingE
540
541 if(secNumberFinal > secNumberInit)
542 {
543 for (std::size_t i=secNumberInit; i<secNumbe
544 {
545 //Check if there is enough residual
546 if (bindingEnergy >= ((*fvect)[i])->
547 {
548 //Ok, this is a valid secondary:
549 bindingEnergy -= ((*fvect)[i])->GetKine
550 //G4cout << "--deex nrj=" << ((*fvect)[
551 //<< G4endl;
552 }
553 else
554 {
555 //Invalid secondary: not enough energy
556 //Keep its energy in the local deposit
557 delete (*fvect)[i];
558 (*fvect)[i]=nullptr;
559 }
560 }
561 }
562
563 //TEST
564 //G4cout << "k=" << k/eV<< G4endl;
565 //G4cout << "secondaryKinetic=" << seconda
566 //G4cout << "scatteredEnergy=" << scattere
567 //G4cout << "deposited energy=" << binding
568 //
569
570 }
571
572 //This should never happen
573 if(bindingEnergy < 0.0)
574 G4Exception("G4DNABornIonisatioModel1::Sa
575 "em2050",FatalException,"Nega
576
577 //bindingEnergy has been decreased
578 //by the amount of energy taken away by de
579 if (!statCode)
580 {
581 fParticleChangeForGamma->SetProposedKine
582 fParticleChangeForGamma->ProposeLocalEne
583 }
584 else
585 {
586 fParticleChangeForGamma->SetProposedKine
587 fParticleChangeForGamma->ProposeLocalEne
588 }
589
590 // TEST //////////////////////////
591 // if (secondaryKinetic<0) abort();
592 // if (scatteredEnergy<0) abort();
593 // if (k-scatteredEnergy-secondaryKinetic-
594 // if (k-scatteredEnergy<0) abort();
595 /////////////////////////////////
596
597 const G4Track * theIncomingTrack = fPartic
598 G4DNAChemistryManager::Instance()->CreateW
599 ionizationShell,
600 theIncomingTrack);
601 }
602 }
603
604 //....oooOO0OOooo........oooOO0OOooo........oo
605
606 G4double G4DNABornIonisationModel1::RandomizeE
607
608
609 {
610 // G4cout << "*** SLOW computation for " <<
611
612 if (particleDefinition == G4Electron::Electr
613 {
614 G4double maximumEnergyTransfer = 0.;
615 if ((k + waterStructure.IonisationEnergy(s
616 maximumEnergyTransfer = k;
617 else
618 maximumEnergyTransfer = (k + waterStruct
619
620 // SI : original method
621 /*
622 G4double crossSectionMaximum = 0.;
623 for(G4double value=waterStructure.Ionisati
624 {
625 G4double differentialCrossSection = Diffe
626 if(differentialCrossSection >= crossSecti
627 }
628 */
629
630 // SI : alternative method
631 G4double crossSectionMaximum = 0.;
632
633 G4double minEnergy = waterStructure.Ionisa
634 G4double maxEnergy = maximumEnergyTransfer
635 G4int nEnergySteps = 50;
636
637 G4double value(minEnergy);
638 G4double stpEnergy(std::pow(maxEnergy / va
639 1. / static_ca
640 G4int step(nEnergySteps);
641 while (step > 0)
642 {
643 step--;
644 G4double differentialCrossSection =
645 DifferentialCrossSection(particleDef
646 k / eV,
647 value / eV,
648 shell);
649 if (differentialCrossSection >= crossSec
650 crossSectionMaximum = differentialCros
651 value *= stpEnergy;
652 }
653 //
654
655 G4double secondaryElectronKineticEnergy =
656 do
657 {
658 secondaryElectronKineticEnergy = G4Unifo
659 } while(G4UniformRand()*crossSectionMaximu
660 DifferentialCrossSection(particleDefin
661 (secondaryElectronKineticEnergy+wa
662
663 return secondaryElectronKineticEnergy;
664
665 }
666
667 if (particleDefinition == G4Proton::ProtonDe
668 {
669 G4double maximumKineticEnergyTransfer = 4.
670 * (electron_mass_c2 / proton_mass_c2)
671
672 G4double crossSectionMaximum = 0.;
673 for (G4double value = waterStructure.Ionis
674 value <= 4. * waterStructure.Ionisatio
675 {
676 G4double differentialCrossSection =
677 DifferentialCrossSection(particleDef
678 k / eV,
679 value / eV,
680 shell);
681 if (differentialCrossSection >= crossSec
682 crossSectionMaximum = differentialCros
683 }
684
685 G4double secondaryElectronKineticEnergy =
686 do
687 {
688 secondaryElectronKineticEnergy = G4Unifo
689 } while(G4UniformRand()*crossSectionMaximu
690 DifferentialCrossSection(particleDefin
691 (secondaryElectronKineticEnergy+wa
692
693 return secondaryElectronKineticEnergy;
694 }
695
696 return 0;
697 }
698
699 //....oooOO0OOooo........oooOO0OOooo........oo
700
701 // The following section is not used anymore b
702 // GetAngularDistribution()->SampleDirectionFo
703
704 /*
705 void G4DNABornIonisationModel1::RandomizeEjec
706 G4double k,
707 G4double secKinetic,
708 G4double & cosTheta,
709 G4double & phi )
710 {
711 if (particleDefinition == G4Electron::Electro
712 {
713 phi = twopi * G4UniformRand();
714 if (secKinetic < 50.*eV) cosTheta = (2.*G4Uni
715 else if (secKinetic <= 200.*eV)
716 {
717 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4
718 else cosTheta = G4UniformRand()*(std::sqrt(2.
719 }
720 else
721 {
722 G4double sin2O = (1.-secKinetic/k) / (1.+secK
723 cosTheta = std::sqrt(1.-sin2O);
724 }
725 }
726
727 else if (particleDefinition == G4Proton::Prot
728 {
729 G4double maxSecKinetic = 4.* (electron_mass_c
730 phi = twopi * G4UniformRand();
731
732 // cosTheta = std::sqrt(secKinetic / maxSecKi
733
734 // Restriction below 100 eV from Emfietzoglou
735
736 if (secKinetic>100*eV) cosTheta = std::sqrt(s
737 else cosTheta = (2.*G4UniformRand())-1.;
738
739 }
740 }
741 */
742
743 //....oooOO0OOooo........oooOO0OOooo........oo
744 G4double G4DNABornIonisationModel1::Differenti
745
746
747
748 {
749 G4double sigma = 0.;
750
751 if (energyTransfer >= waterStructure.Ionisat
752 {
753 G4double valueT1 = 0;
754 G4double valueT2 = 0;
755 G4double valueE21 = 0;
756 G4double valueE22 = 0;
757 G4double valueE12 = 0;
758 G4double valueE11 = 0;
759
760 G4double xs11 = 0;
761 G4double xs12 = 0;
762 G4double xs21 = 0;
763 G4double xs22 = 0;
764
765 if (particleDefinition == G4Electron::Elec
766 {
767
768 // Protection against out of boundary ac
769 if (k==eTdummyVec.back()) k=k*(1.-1e-12)
770 //
771
772 // k should be in eV and energy transfer
773
774 auto t2 = std::upper_bound(eTdummyVec.be
775
776
777
778 auto t1 = t2 - 1;
779
780 // SI : the following condition avoids s
781 if (energyTransfer <= eVecm[(*t1)].back(
782 && energyTransfer <= eVecm[(*t2)].ba
783 {
784 auto e12 =
785 std::upper_bound(eVecm[(*t1)].begi
786 eVecm[(*t1)].end(
787 energyTransfer);
788 auto e11 = e12 - 1;
789
790 auto e22 =
791 std::upper_bound(eVecm[(*t2)].begi
792 eVecm[(*t2)].end(
793 energyTransfer);
794 auto e21 = e22 - 1;
795
796 valueT1 = *t1;
797 valueT2 = *t2;
798 valueE21 = *e21;
799 valueE22 = *e22;
800 valueE12 = *e12;
801 valueE11 = *e11;
802
803 xs11 = eDiffCrossSectionData[ionizatio
804 xs12 = eDiffCrossSectionData[ionizatio
805 xs21 = eDiffCrossSectionData[ionizatio
806 xs22 = eDiffCrossSectionData[ionizatio
807
808 }
809
810 }
811
812 if (particleDefinition == G4Proton::Proton
813 {
814 // Protection against out of boundary ac
815 if (k==pTdummyVec.back()) k=k*(1.-1e-12)
816 //
817
818 // k should be in eV and energy transfer
819 auto t2 = std::upper_bound(pTdummyVec.be
820
821
822 auto t1 = t2 - 1;
823
824 auto e12 = std::upper_bound(pVecm[(*t1)]
825
826
827 auto e11 = e12 - 1;
828
829 auto e22 = std::upper_bound(pVecm[(*t2)]
830
831
832 auto e21 = e22 - 1;
833
834 valueT1 = *t1;
835 valueT2 = *t2;
836 valueE21 = *e21;
837 valueE22 = *e22;
838 valueE12 = *e12;
839 valueE11 = *e11;
840
841 xs11 = pDiffCrossSectionData[ionizationL
842 xs12 = pDiffCrossSectionData[ionizationL
843 xs21 = pDiffCrossSectionData[ionizationL
844 xs22 = pDiffCrossSectionData[ionizationL
845
846 }
847
848 G4double xsProduct = xs11 * xs12 * xs21 *
849 if (xsProduct != 0.)
850 {
851 sigma = QuadInterpolator(valueE11,
852 valueE12,
853 valueE21,
854 valueE22,
855 xs11,
856 xs12,
857 xs21,
858 xs22,
859 valueT1,
860 valueT2,
861 k,
862 energyTransfer)
863 }
864 }
865
866 return sigma;
867 }
868
869 //....oooOO0OOooo........oooOO0OOooo........oo
870
871 G4double G4DNABornIonisationModel1::Interpolat
872
873
874
875
876 {
877 G4double value = 0.;
878
879 // Log-log interpolation by default
880
881 if (e1 != 0 && e2 != 0 && (std::log10(e2) -
882 && !fasterCode)
883 {
884 G4double a = (std::log10(xs2) - std::log10
885 / (std::log10(e2) - std::log10(e1));
886 G4double b = std::log10(xs2) - a * std::lo
887 G4double sigma = a * std::log10(e) + b;
888 value = (std::pow(10., sigma));
889 }
890
891 // Switch to lin-lin interpolation
892 /*
893 if ((e2-e1)!=0)
894 {
895 G4double d1 = xs1;
896 G4double d2 = xs2;
897 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1)
898 }
899 */
900
901 // Switch to log-lin interpolation for faste
902 if ((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 &
903 {
904 G4double d1 = std::log10(xs1);
905 G4double d2 = std::log10(xs2);
906 value = std::pow(10., (d1 + (d2 - d1) * (e
907 }
908
909 // Switch to lin-lin interpolation for faste
910 // in case one of xs1 or xs2 (=cum proba) va
911
912 if ((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0)
913 {
914 G4double d1 = xs1;
915 G4double d2 = xs2;
916 value = (d1 + (d2 - d1) * (e - e1) / (e2 -
917 }
918
919 /*
920 G4cout
921 << e1 << " "
922 << e2 << " "
923 << e << " "
924 << xs1 << " "
925 << xs2 << " "
926 << value
927 << G4endl;
928 */
929
930 return value;
931 }
932
933 //....oooOO0OOooo........oooOO0OOooo........oo
934
935 G4double G4DNABornIonisationModel1::QuadInterp
936
937
938
939
940
941
942
943
944
945
946
947 {
948 G4double interpolatedvalue1 = Interpolate(e1
949 G4double interpolatedvalue2 = Interpolate(e2
950 G4double value = Interpolate(t1,
951 t2,
952 t,
953 interpolatedval
954 interpolatedval
955
956 return value;
957 }
958
959 G4double G4DNABornIonisationModel1::GetPartial
960
961
962
963 {
964 std::map<G4String, G4DNACrossSectionDataSet*
965 pos = tableData.find(particle->GetParticleNa
966
967 if (pos != tableData.end())
968 {
969 G4DNACrossSectionDataSet* table = pos->sec
970 return table->GetComponent(level)->FindVal
971 }
972
973 return 0;
974 }
975
976 //....oooOO0OOooo........oooOO0OOooo........oo
977
978 G4int G4DNABornIonisationModel1::RandomSelect(
979
980 {
981 G4int level = 0;
982
983 std::map<G4String, G4DNACrossSectionDataSet*
984 pos = tableData.find(particle);
985
986 if (pos != tableData.end())
987 {
988 G4DNACrossSectionDataSet* table = pos->sec
989
990 if (table != nullptr)
991 {
992 auto valuesBuffer = new G4double[table-
993 const auto n = (G4int)table->NumberOfCo
994 G4int i(n);
995 G4double value = 0.;
996
997 while (i > 0)
998 {
999 i--;
1000 valuesBuffer[i] = table->GetComponent
1001 value += valuesBuffer[i];
1002 }
1003
1004 value *= G4UniformRand();
1005
1006 i = n;
1007
1008 while (i > 0)
1009 {
1010 i--;
1011
1012 if (valuesBuffer[i] > value)
1013 {
1014 delete[] valuesBuffer;
1015 return i;
1016 }
1017 value -= valuesBuffer[i];
1018 }
1019
1020
1021 delete[] valuesBuffer;
1022
1023 }
1024 } else
1025 {
1026 G4Exception("G4DNABornIonisationModel1::R
1027 "em0002",
1028 FatalException,
1029 "Model not applicable to part
1030 }
1031
1032 return level;
1033 }
1034
1035 //....oooOO0OOooo........oooOO0OOooo........o
1036
1037 G4double G4DNABornIonisationModel1::Randomize
1038
1039
1040 {
1041 //G4cout << "*** FAST computation for " <<
1042
1043 G4double secondaryElectronKineticEnergy = 0
1044
1045 G4double random = G4UniformRand();
1046
1047 secondaryElectronKineticEnergy = Transfered
1048
1049
1050
1051 - waterStructure.IonisationEnergy(shell
1052
1053 //G4cout << RandomTransferedEnergy(particle
1054
1055 if (secondaryElectronKineticEnergy < 0.)
1056 return 0.;
1057 //
1058
1059 return secondaryElectronKineticEnergy;
1060 }
1061
1062 //....oooOO0OOooo........oooOO0OOooo........o
1063
1064 G4double G4DNABornIonisationModel1::Transfere
1065
1066
1067
1068 {
1069 G4double nrj = 0.;
1070
1071 G4double valueK1 = 0;
1072 G4double valueK2 = 0;
1073 G4double valuePROB21 = 0;
1074 G4double valuePROB22 = 0;
1075 G4double valuePROB12 = 0;
1076 G4double valuePROB11 = 0;
1077
1078 G4double nrjTransf11 = 0;
1079 G4double nrjTransf12 = 0;
1080 G4double nrjTransf21 = 0;
1081 G4double nrjTransf22 = 0;
1082
1083 if (particleDefinition == G4Electron::Elect
1084 {
1085 // Protection against out of boundary acc
1086 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
1087 //
1088
1089 // k should be in eV
1090 auto k2 = std::upper_bound(eTdummyVec.beg
1091
1092
1093 auto k1 = k2 - 1;
1094
1095 /*
1096 G4cout << "----> k=" << k
1097 << " " << *k1
1098 << " " << *k2
1099 << " " << random
1100 << " " << ionizationLevelIndex
1101 << " " << eProbaShellMap[ionizationLevel
1102 << " " << eProbaShellMap[ionizationLevel
1103 << G4endl;
1104 */
1105
1106 // SI : the following condition avoids si
1107 if (random <= eProbaShellMap[ionizationLe
1108 && random <= eProbaShellMap[ionizatio
1109 {
1110 auto prob12 =
1111 std::upper_bound(eProbaShellMap[ion
1112 eProbaShellMap[ion
1113 random);
1114
1115 auto prob11 = prob12 - 1;
1116
1117 auto prob22 =
1118 std::upper_bound(eProbaShellMap[ion
1119 eProbaShellMap[ion
1120 random);
1121
1122 auto prob21 = prob22 - 1;
1123
1124 valueK1 = *k1;
1125 valueK2 = *k2;
1126 valuePROB21 = *prob21;
1127 valuePROB22 = *prob22;
1128 valuePROB12 = *prob12;
1129 valuePROB11 = *prob11;
1130
1131 /*
1132 G4cout << " " << random << " "
1133 << valuePROB12 << " " << valuePROB21 <
1134 */
1135
1136 nrjTransf11 = eNrjTransfData[ionization
1137 nrjTransf12 = eNrjTransfData[ionization
1138 nrjTransf21 = eNrjTransfData[ionization
1139 nrjTransf22 = eNrjTransfData[ionization
1140
1141 /*
1142 G4cout << " " << ionizationLeve
1143 << random << " " <<valueK1 << " " << v
1144
1145 G4cout << " " << random << " "
1146 << nrjTransf12 << " " << nrjTransf21 <
1147 */
1148
1149 }
1150
1151 // Avoids cases where cum xs is zero for
1152 if (random > eProbaShellMap[ionizationLev
1153 {
1154 auto prob22 =
1155 std::upper_bound(eProbaShellMap[ion
1156 eProbaShellMap[ion
1157 random);
1158
1159 auto prob21 = prob22 - 1;
1160
1161 valueK1 = *k1;
1162 valueK2 = *k2;
1163 valuePROB21 = *prob21;
1164 valuePROB22 = *prob22;
1165
1166 //G4cout << " " << random << " "
1167
1168 nrjTransf21 = eNrjTransfData[ionization
1169 nrjTransf22 = eNrjTransfData[ionization
1170
1171 G4double interpolatedvalue2 = Interpola
1172
1173
1174
1175
1176
1177 // zeros are explicitly set
1178
1179 G4double value = Interpolate(valueK1, v
1180
1181 /*
1182 G4cout << " " << ionizationLeve
1183 << random << " " <<valueK1 << " " << v
1184
1185 G4cout << " " << random << " "
1186 << nrjTransf12 << " " << nrjTransf21 <
1187
1188 G4cout << "ici" << " " << value << G4e
1189 */
1190
1191 return value;
1192 }
1193 }
1194 //
1195 else if (particleDefinition == G4Proton::Pr
1196 {
1197 // Protection against out of boundary acc
1198 if (k==pTdummyVec.back()) k=k*(1.-1e-12);
1199 //
1200
1201 // k should be in eV
1202
1203 auto k2 = std::upper_bound(pTdummyVec.beg
1204
1205
1206
1207 auto k1 = k2 - 1;
1208
1209 /*
1210 G4cout << "----> k=" << k
1211 << " " << *k1
1212 << " " << *k2
1213 << " " << random
1214 << " " << ionizationLevelIndex
1215 << " " << pProbaShellMap[ionizationLevel
1216 << " " << pProbaShellMap[ionizationLevel
1217 << G4endl;
1218 */
1219
1220 // SI : the following condition avoids si
1221 // for eg. when the last element is
1222 if (random <= pProbaShellMap[ionizationLe
1223 && random <= pProbaShellMap[ionizatio
1224 {
1225 auto prob12 =
1226 std::upper_bound(pProbaShellMap[ion
1227 pProbaShellMap[ion
1228 random);
1229
1230 auto prob11 = prob12 - 1;
1231
1232 auto prob22 =
1233 std::upper_bound(pProbaShellMap[ion
1234 pProbaShellMap[ion
1235 random);
1236
1237 auto prob21 = prob22 - 1;
1238
1239 valueK1 = *k1;
1240 valueK2 = *k2;
1241 valuePROB21 = *prob21;
1242 valuePROB22 = *prob22;
1243 valuePROB12 = *prob12;
1244 valuePROB11 = *prob11;
1245
1246 /*
1247 G4cout << " " << random << " "
1248 << valuePROB12 << " " << valuePROB21 <
1249 */
1250
1251 nrjTransf11 = pNrjTransfData[ionization
1252 nrjTransf12 = pNrjTransfData[ionization
1253 nrjTransf21 = pNrjTransfData[ionization
1254 nrjTransf22 = pNrjTransfData[ionization
1255
1256 /*
1257 G4cout << " " << ionizationLeve
1258 << random << " " <<valueK1 << " " << v
1259
1260 G4cout << " " << random << " "
1261 << nrjTransf12 << " " << nrjTransf21 <
1262 */
1263 }
1264
1265 // Avoids cases where cum xs is zero for
1266
1267 if (random > pProbaShellMap[ionizationLev
1268 {
1269 auto prob22 =
1270 std::upper_bound(pProbaShellMap[ion
1271 pProbaShellMap[ion
1272 random);
1273
1274 auto prob21 = prob22 - 1;
1275
1276 valueK1 = *k1;
1277 valueK2 = *k2;
1278 valuePROB21 = *prob21;
1279 valuePROB22 = *prob22;
1280
1281 //G4cout << " " << random << " "
1282
1283 nrjTransf21 = pNrjTransfData[ionization
1284 nrjTransf22 = pNrjTransfData[ionization
1285
1286 G4double interpolatedvalue2 = Interpola
1287
1288
1289
1290
1291
1292 // zeros are explicitly set
1293
1294 G4double value = Interpolate(valueK1, v
1295
1296 /*
1297 G4cout << " " << ionizationLeve
1298 << random << " " <<valueK1 << " " << v
1299
1300 G4cout << " " << random << " "
1301 << nrjTransf12 << " " << nrjTransf21 <
1302
1303 G4cout << "ici" << " " << value << G4e
1304 */
1305
1306 return value;
1307 }
1308 }
1309 // End electron and proton cases
1310
1311 G4double nrjTransfProduct = nrjTransf11 * n
1312 * nrjTransf22;
1313 //G4cout << "nrjTransfProduct=" << nrjTrans
1314
1315 if (nrjTransfProduct != 0.)
1316 {
1317 nrj = QuadInterpolator(valuePROB11,
1318 valuePROB12,
1319 valuePROB21,
1320 valuePROB22,
1321 nrjTransf11,
1322 nrjTransf12,
1323 nrjTransf21,
1324 nrjTransf22,
1325 valueK1,
1326 valueK2,
1327 k,
1328 random);
1329 }
1330 //G4cout << nrj << endl;
1331
1332 return nrj;
1333 }
1334