| 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 // $Id: G4DNARelativisticIonisationModel.cc $
27 //
28 // Created on 2016/05/12
29 //
30 // Authors: D Sakata, S. Incerti
31 //
32 // This class perform ionisation for electron transportation in gold,
33 // based on Relativistic Binary Encounter Bethe-Vriens(RBEBV) model.
34 // See following reference paper,
35 // M. Guerra et al, J. Phys. B: At. Mol. Opt. Phys. 48, 185202 (2015)
36 // =======================================================================
37 // Limitation of secondaries by GEANT4 atomic de-excitation:
38 // The cross section and energy of secondary production is based on
39 // EADL database. If there are no tabele for several orbitals, this class
40 // will not provide secondaries for the orbitals.
41 // For gold(Au), this class provide secondaries for inner 18 orbitals
42 // but don't provide for outer 3 orbitals due to EADL databese limitation.
43 // =======================================================================
44
45 #include "G4DNARelativisticIonisationModel.hh"
46 #include "G4SystemOfUnits.hh"
47 #include "G4AtomicShell.hh"
48 #include "G4UAtomicDeexcitation.hh"
49 #include "G4LossTableManager.hh"
50 #include "G4Gamma.hh"
51 #include "G4RandomDirection.hh"
52
53 #include "G4DNAMolecularMaterial.hh"
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56
57 using namespace std;
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61 G4DNARelativisticIonisationModel::G4DNARelativisticIonisationModel(
62 const G4ParticleDefinition*,
63 const G4String& nam) :
64 G4VEmModel(nam), fasterCode(true)
65 {
66 fHighEnergyLimit = 0;
67 fLowEnergyLimit = 0;
68
69 verboseLevel = 0;
70
71 SetDeexcitationFlag(true);
72 fAtomDeexcitation = nullptr;
73 fMaterialDensity = nullptr;
74 fParticleDefinition = nullptr;
75 fParticleChangeForGamma = nullptr;
76
77 if (verboseLevel > 0)
78 {
79 G4cout << "Relativistic Ionisation Model is constructed " << G4endl;
80 }
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85 G4DNARelativisticIonisationModel::~G4DNARelativisticIonisationModel()
86 = default;
87
88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
89
90 void G4DNARelativisticIonisationModel::Initialise(const G4ParticleDefinition* particle,
91 const G4DataVector& /*cuts*/)
92 {
93
94 if (verboseLevel > 3)
95 {
96 G4cout <<
97 "Calling G4DNARelativisticIonisationModel::Initialise()"
98 << G4endl;
99 }
100
101
102 if(fParticleDefinition != nullptr && fParticleDefinition != particle)
103 {
104 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0001",
105 FatalException,"Model already initialized for another particle type.");
106 }
107
108 fParticleDefinition = particle;
109 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
110 if(particle == electronDef)
111 {
112 fLowEnergyLimit = 10 * eV;
113 fHighEnergyLimit = 1.0 * GeV;
114
115 std::ostringstream eFullFileNameZ;
116
117 const char *path = G4FindDataDir("G4LEDATA");
118 if (path == nullptr)
119 {
120 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0006",
121 FatalException,"G4LEDATA environment variable not set.");
122 return;
123 }
124
125
126 G4ProductionCutsTable *coupletable
127 = G4ProductionCutsTable::GetProductionCutsTable();
128 auto Ncouple = (G4int)coupletable ->GetTableSize();
129 for(G4int i=0; i<Ncouple; ++i)
130 {
131 const G4MaterialCutsCouple* couple
132 = coupletable->GetMaterialCutsCouple(i);
133 const G4Material * material = couple ->GetMaterial();
134 {
135 // Protection: only for single element
136 if(material->GetNumberOfElements()>1) continue;
137
138 G4int Z = material->GetZ();
139 // Protection: only for GOLD
140 if(Z!=79) continue;
141
142 iState [Z].clear();
143 iShell [Z].clear();
144 iSubShell [Z].clear();
145 Nelectrons[Z].clear();
146 Ebinding [Z].clear();
147 Ekinetic [Z].clear();
148 LoadAtomicStates(Z,path);
149
150 /////////////Load cumulated DCS////////////////
151 eVecEZ.clear();
152 eVecEjeEZ.clear();
153 eProbaShellMapZ.clear();
154 eDiffCrossSectionDataZ.clear();
155
156 eFullFileNameZ.str("");
157 eFullFileNameZ.clear(stringstream::goodbit);
158
159 eFullFileNameZ
160 << path
161 << "/dna/sigmadiff_cumulated_ionisation_e_RBEBV_Z"
162 << Z << ".dat";
163 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
164 if (!eDiffCrossSectionZ)
165 G4Exception("G4DNARelativisticIonisationModel::Initialise","em0003",
166 FatalException,
167 "Missing data file for cumulated DCS");
168
169 eVecEZ[Z].push_back(0.);
170 while(!eDiffCrossSectionZ.eof())
171 {
172 G4double tDummy;
173 G4double eDummy;
174 eDiffCrossSectionZ>>tDummy>>eDummy;
175 if (tDummy != eVecEZ[Z].back())
176 {
177 eVecEZ[Z].push_back(tDummy);
178 eVecEjeEZ[Z][tDummy].push_back(0.);
179 }
180
181 for(G4int istate=0;istate<(G4int)iState[Z].size();istate++)
182 {
183 eDiffCrossSectionZ>>
184 eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy];
185 eEjectedEnergyDataZ[Z][istate][tDummy]
186 [eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy]]
187 = eDummy;
188 eProbaShellMapZ[Z][istate][tDummy].push_back(
189 eDiffCrossSectionDataZ[Z][istate][tDummy][eDummy]);
190 }
191
192 if (eDummy != eVecEjeEZ[Z][tDummy].back()){
193 eVecEjeEZ[Z][tDummy].push_back(eDummy);
194 }
195 }
196 }
197 }
198 }
199 else
200 {
201 G4cout<<
202 "Error : No particle Definition is found in G4DNARelativisticIonisationModel"
203 <<G4endl;
204 return;
205 }
206
207 if( verboseLevel>0 )
208 {
209 G4cout << "Relativistic Ionisation model is initialized " << G4endl
210 << "Energy range: "
211 << LowEnergyLimit() / eV << " eV - "
212 << HighEnergyLimit() / keV << " keV for "
213 << particle->GetParticleName()
214 << G4endl;
215 }
216
217 // Initialise gold density pointer
218 fMaterialDensity = G4DNAMolecularMaterial::Instance()
219 ->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_Au"));
220
221 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
222 fParticleChangeForGamma = GetParticleChangeForGamma();
223
224 if (isInitialised){return;}
225 isInitialised = true;
226 }
227
228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
229
230 G4double G4DNARelativisticIonisationModel::CrossSectionPerVolume(
231 const G4Material* material,
232 const G4ParticleDefinition* particleDefinition,
233 G4double ekin,
234 G4double,
235 G4double)
236 {
237 if (verboseLevel > 3)
238 {
239 G4cout <<
240 "Calling CrossSectionPerVolume() of G4DNARelativisticIonisationModel"
241 << G4endl;
242 }
243
244 if(particleDefinition != fParticleDefinition) return 0;
245
246 // Calculate total cross section for model
247 G4double sigma=0;
248
249 if(material->GetNumberOfElements()>1) return 0.; // Protection for Molecules
250 G4double atomicNDensity = material->GetAtomicNumDensityVector()[0];
251 G4double z = material->GetZ();
252
253 if(atomicNDensity!= 0.0)
254 {
255 if (ekin >= fLowEnergyLimit && ekin < fHighEnergyLimit)
256 {
257 sigma = GetTotalCrossSection(material,particleDefinition,ekin);
258 }
259
260 if (verboseLevel > 2)
261 {
262 G4cout << "__________________________________" << G4endl;
263 G4cout << "=== G4DNARelativisticIonisationModel - XS INFO START" <<G4endl;
264 G4cout << "=== Kinetic energy (eV)=" << ekin/eV << " particle : "
265 << particleDefinition->GetParticleName() << G4endl;
266 G4cout << "=== Cross section per atom for Z="<<z<<" is (cm^2)"
267 << sigma/cm/cm << G4endl;
268 G4cout << "=== Cross section per atom for Z="<<z<<" is (cm^-1)="
269 << sigma*atomicNDensity/(1./cm) << G4endl;
270 G4cout << "=== G4DNARelativisticIonisationModel - XS INFO END" << G4endl;
271 }
272 }
273 return sigma*atomicNDensity;
274 }
275
276 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277
278 void G4DNARelativisticIonisationModel::SampleSecondaries(
279 std::vector<G4DynamicParticle*>* fvect,
280 const G4MaterialCutsCouple* couple,
281 const G4DynamicParticle* particle,
282 G4double,G4double)
283 {
284 if (verboseLevel > 3)
285 {
286 G4cout <<
287 "Calling SampleSecondaries() of G4DNARelativisticIonisationModel"
288 << G4endl;
289 }
290
291
292 G4ParticleDefinition* particleDef = particle->GetDefinition();
293 G4double k = particle->GetKineticEnergy();
294 G4double ejectedE = 0.*eV;
295
296 if(fLowEnergyLimit <= k && k<fHighEnergyLimit)
297 {
298 G4ThreeVector primaryDir = particle ->GetMomentumDirection();
299
300 G4double particleMass = particleDef->GetPDGMass();
301 G4double totalEnergy = k+particleMass;
302 G4double pSquare = k*(totalEnergy+particleMass);
303 G4double totalMomentum = std::sqrt(pSquare);
304
305 const G4Material *material = couple->GetMaterial();
306 G4int z = material->GetZ();
307 G4int level = RandomSelect(material,particleDef,k);
308
309 if(k<Ebinding[z].at(level)) return;
310
311 G4int NumSecParticlesInit =0;
312 G4int NumSecParticlesFinal=0;
313
314 if(fAtomDeexcitation != nullptr){
315 auto as = G4AtomicShellEnumerator(level);
316 const G4AtomicShell *shell = fAtomDeexcitation->GetAtomicShell(z,as);
317 NumSecParticlesInit = (G4int)fvect->size();
318 fAtomDeexcitation->GenerateParticles(fvect,shell,z,0,0);
319 NumSecParticlesFinal = (G4int)fvect->size();
320 }
321
322 ejectedE
323 = GetEjectedElectronEnergy (material,particleDef,k,level);
324 G4ThreeVector ejectedDir
325 = GetEjectedElectronDirection(particleDef,k,ejectedE);
326 ejectedDir.rotateUz(primaryDir);
327
328 G4double scatteredE = k - Ebinding[z].at(level) - ejectedE;
329
330 if(particleDef == G4Electron::ElectronDefinition()){
331 G4double secondaryTotMomentum
332 = std::sqrt(ejectedE*(ejectedE+2*CLHEP::electron_mass_c2));
333 G4double finalMomentumX
334 = totalMomentum*primaryDir.x()- secondaryTotMomentum*ejectedDir.x();
335 G4double finalMomentumY
336 = totalMomentum*primaryDir.y()- secondaryTotMomentum*ejectedDir.y();
337 G4double finalMomentumZ
338 = totalMomentum*primaryDir.z()- secondaryTotMomentum*ejectedDir.z();
339
340 G4ThreeVector scatteredDir(finalMomentumX,finalMomentumY,finalMomentumZ);
341 fParticleChangeForGamma->ProposeMomentumDirection(scatteredDir.unit());
342
343 }
344 else
345 {
346 fParticleChangeForGamma->ProposeMomentumDirection(primaryDir);
347 }
348
349 //G4double deexSecEnergy=0.;
350 G4double restEproduction = Ebinding[z].at(level);
351 for(G4int iparticle=NumSecParticlesInit;
352 iparticle<NumSecParticlesFinal;iparticle++)
353 {
354 //deexSecEnergy = deexSecEnergy + (*fvect)[iparticle]->GetKineticEnergy();
355 G4double Edeex = (*fvect)[iparticle]->GetKineticEnergy();
356 if(restEproduction>=Edeex){
357 restEproduction -= Edeex;
358 }
359 else{
360 delete (*fvect)[iparticle];
361 (*fvect)[iparticle]=nullptr;
362 }
363 }
364 if(restEproduction < 0.0){
365 G4Exception("G4DNARelativisticIonisationModel::SampleSecondaries()",
366 "em0008",FatalException,"Negative local energy deposit");
367 }
368
369 if(!statCode)
370 {
371 if(scatteredE>0){
372 fParticleChangeForGamma->SetProposedKineticEnergy (scatteredE);
373 fParticleChangeForGamma->ProposeLocalEnergyDeposit(restEproduction);
374 //fParticleChangeForGamma
375 //->ProposeLocalEnergyDeposit(k-scatteredE-ejectedE-deexSecEnergy);
376 }
377 }
378 else
379 {
380 fParticleChangeForGamma->SetProposedKineticEnergy (k);
381 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredE);
382 }
383
384 if(ejectedE>0){
385 auto ejectedelectron
386 = new G4DynamicParticle(G4Electron::Electron(),ejectedDir,ejectedE);
387 fvect->push_back(ejectedelectron);
388 }
389 }
390 }
391
392 void G4DNARelativisticIonisationModel::LoadAtomicStates(
393 G4int z,const char* path)
394 {
395
396 if (verboseLevel > 3)
397 {
398 G4cout <<
399 "Calling LoadAtomicStates() of G4DNARelativisticIonisationModel"
400 << G4endl;
401 }
402 const char *datadir = path;
403 if(datadir == nullptr)
404 {
405 datadir = G4FindDataDir("G4LEDATA");
406 if(datadir == nullptr)
407 {
408 G4Exception("G4DNARelativisticIonisationModel::LoadAtomicStates()",
409 "em0002",FatalException,"Enviroment variable G4LEDATA not defined");
410
411 return;
412 }
413 }
414 std::ostringstream targetfile;
415 targetfile << datadir <<"/dna/atomicstate_Z"<< z <<".dat";
416 std::ifstream fin(targetfile.str().c_str());
417 if(!fin)
418 {
419 G4cout<< " Error : "<< targetfile.str() <<" is not found "<<G4endl;
420 G4Exception("G4DNARelativisticIonisationModel::LoadAtomicStates()","em0002",
421 FatalException,"There is no target file");
422 return;
423 }
424
425 G4String buff0,buff1,buff2,buff3,buff4,buff5,buff6;
426 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
427 G4int iline=0;
428 while(true){
429 fin >> buff0 >>buff1>>buff2>>buff3>>buff4>>buff5>>buff6;
430 if(!fin.eof())
431 {
432 iState [z].push_back(stoi(buff0));
433 iShell [z].push_back(stoi(buff1));
434 iSubShell [z].push_back(stoi(buff2));
435 Nelectrons[z].push_back(stoi(buff3));
436 Ebinding [z].push_back(stod(buff4));
437 if(stod(buff5)==0.)
438 {// if there is no kinetic energy in the file, kinetic energy
439 // for Bhor atomic model will be calculated: !!! I's not realistic!!!
440 G4double radius = std::pow(iShell[z].at(iline),2)
441 *std::pow(CLHEP::hbar_Planck,2)*(4*CLHEP::pi*CLHEP::epsilon0)
442 /CLHEP::electron_mass_c2;
443 G4double momentum = iShell[z].at(iline)*CLHEP::hbar_Planck/radius;
444 Ekinetic[z].push_back(std::pow(momentum,2)/(2*CLHEP::electron_mass_c2));
445 }
446 else
447 {
448 Ekinetic [z].push_back(stod(buff5));
449 }
450 iline++;
451 }
452 else
453 {
454 break;
455 }
456 }
457 }
458
459 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
460
461 G4double G4DNARelativisticIonisationModel::GetTotalCrossSection(
462 const G4Material* material,
463 const G4ParticleDefinition* particle,
464 G4double kineticEnergy)
465 {
466 G4double value=0;
467 G4int z = material->GetZ();
468 if(z!=79){ return 0.;}
469
470 std::size_t N=iState[z].size();
471 for(G4int i=0; i<(G4int)N; ++i){
472 value = value+GetPartialCrossSection(material,i,particle,kineticEnergy);
473 }
474 return value;
475 }
476
477 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
478
479 G4double G4DNARelativisticIonisationModel::GetPartialCrossSection(
480 const G4Material* material,
481 G4int level,
482 const G4ParticleDefinition* particle,
483 G4double kineticEnergy)
484 {
485 G4double value = 0;
486 G4double constRy =13.6057E-6;//MeV
487
488 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
489 G4int z = material->GetZ();
490 if(particle==electronDef){
491
492 G4double t = kineticEnergy /Ebinding[z].at(level);
493 G4double tdash = kineticEnergy /CLHEP::electron_mass_c2;
494 G4double udash = Ekinetic[z].at(level)/CLHEP::electron_mass_c2;
495 G4double bdash = Ebinding[z].at(level)/CLHEP::electron_mass_c2;
496 G4double beta_t2 = 1.-1./std::pow(1.+tdash,2);
497 G4double beta_u2 = 1.-1./std::pow(1.+udash,2);
498 G4double beta_b2 = 1.-1./std::pow(1.+bdash,2);
499 G4double alpha = std::sqrt(2*constRy/CLHEP::electron_mass_c2);
500 G4double phi = std::cos(std::sqrt(std::pow(alpha,2)
501 /(beta_t2+beta_b2))*G4Log(beta_t2/beta_b2));
502 G4double constS = 4*CLHEP::pi*std::pow(CLHEP::Bohr_radius,2)
503 *Nelectrons[z].at(level)*std::pow(alpha,4);
504
505 if(Ebinding[z].at(level)<=kineticEnergy)
506 {
507 value =constS/((beta_t2+(beta_u2+beta_b2)/iShell[z].at(level))*2.*bdash)
508 *(1./2.*(G4Log(beta_t2/(1.-beta_t2))-beta_t2-G4Log(2.*bdash))
509 *(1.-1./std::pow(t,2.))
510 +1.-1./t-G4Log(t)/(t+1.)*(1.+2.*tdash)/(std::pow(1.+tdash/2.,2.))
511 *phi+std::pow(bdash,2)/(std::pow(1+tdash/2.,2))*(t-1)/2.);
512 }
513
514 }
515 return value;
516 }
517
518 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
519
520 G4double G4DNARelativisticIonisationModel::GetDifferentialCrossSection(
521 const G4Material* material,
522 const G4ParticleDefinition* particle,
523 G4double kineticEnergy,
524 G4double secondaryEnergy,
525 G4int level)
526 {
527 G4double value=0.;
528 G4double constRy =13.6057E-6;//MeV
529
530 G4int z = material->GetZ();
531
532 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
533 if(particle==electronDef){
534 G4double w = secondaryEnergy /Ebinding[z].at(level);
535 G4double t = kineticEnergy /Ebinding[z].at(level);
536 G4double tdash = kineticEnergy /CLHEP::electron_mass_c2;
537 G4double udash = Ekinetic[z].at(level)/CLHEP::electron_mass_c2;
538 G4double bdash = Ebinding[z].at(level)/CLHEP::electron_mass_c2;
539 G4double beta_t2 = 1.-1./std::pow(1.+tdash,2);
540 G4double beta_u2 = 1.-1./std::pow(1.+udash,2);
541 G4double beta_b2 = 1.-1./std::pow(1.+bdash,2);
542 G4double alpha = std::sqrt(2*constRy/CLHEP::electron_mass_c2);
543 G4double phi = std::cos(std::sqrt(std::pow(alpha,2)/(beta_t2+beta_b2))
544 *G4Log(beta_t2/beta_b2));
545 G4double constS = 4*CLHEP::pi*std::pow(CLHEP::Bohr_radius,2)
546 *Nelectrons[z].at(level)*std::pow(alpha,4);
547
548 if(secondaryEnergy<=((kineticEnergy-Ebinding[z].at(level))/2.))
549 {
550 value = constS/((beta_t2+(beta_u2+beta_b2)/iShell[z].at(level))*2.*bdash)
551 *(-phi/(t+1.)*(1./std::pow(w+1.,1.)+1./std::pow(t-w,1.))
552 *(1.+2*tdash)/std::pow(1.+tdash/2.,2.)
553 +1./std::pow(w+1.,2.)+1./std::pow(t-w,2.)
554 +std::pow(bdash,2)/std::pow(1+tdash/2.,2)
555 +(1./std::pow(w+1.,3.)+1./std::pow(t-w,3.))
556 *(G4Log(beta_t2/(1.-beta_t2))-beta_t2-G4Log(2*bdash)));
557 }
558 }
559 return value;
560 }
561
562 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
563
564 G4int G4DNARelativisticIonisationModel::RandomSelect(
565 const G4Material* material,
566 const G4ParticleDefinition* particle,
567 G4double kineticEnergy)
568 {
569 G4double value = 0.;
570 G4int z = material->GetZ();
571 std::size_t numberOfShell = iShell[z].size();
572 std::vector<G4double> valuesBuffer(numberOfShell, 0.0);
573 const auto n = (G4int)iShell[z].size();
574 G4int i(n);
575
576 while (i > 0)
577 {
578 i--;
579 if((fLowEnergyLimit<=kineticEnergy)&&(kineticEnergy<fHighEnergyLimit))
580 {
581 valuesBuffer[i]=GetPartialCrossSection(material,i,particle,kineticEnergy);
582 }
583 value += valuesBuffer[i];
584 }
585
586 value *= G4UniformRand();
587 i = n;
588
589 while (i > 0)
590 {
591 i--;
592
593 if (valuesBuffer[i] > value)
594 {
595 return i;
596 }
597 value -= valuesBuffer[i];
598 }
599
600 return 9999;
601 }
602
603 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
604
605
606 G4double G4DNARelativisticIonisationModel::GetEjectedElectronEnergy(
607 const G4Material* material,
608 const G4ParticleDefinition* particle,
609 G4double energy, G4int ishell)
610 {
611 G4double secondaryEnergy=0;
612
613 G4ParticleDefinition *electronDef = G4Electron::ElectronDefinition();
614 G4int z = material->GetZ();
615 if(!fasterCode){ // for 2D rejection method
616 if(particle==electronDef){
617 G4double maximumsecondaryEnergy = (energy-Ebinding[z].at(ishell))/2.;
618 if(maximumsecondaryEnergy<0.) return 0.;
619 G4double maximumCrossSection=-999.;
620
621 maximumCrossSection
622 = GetDifferentialCrossSection(material,particle,energy,0.,ishell);
623 do{
624 secondaryEnergy = G4UniformRand()* maximumsecondaryEnergy;
625 }while(G4UniformRand()*maximumCrossSection >
626 GetDifferentialCrossSection(
627 material,particle,energy,secondaryEnergy,ishell));
628 }
629 }
630 else { // for cumulative method using cumulated DCS file
631
632 G4double valueE1 =0.;
633 G4double valueE2 =0.;
634 G4double valueXS21=0.;
635 G4double valueXS22=0.;
636 G4double valueXS11=0.;
637 G4double valueXS12=0.;
638 G4double ejeE21 =0.;
639 G4double ejeE22 =0.;
640 G4double ejeE11 =0.;
641 G4double ejeE12 =0.;
642 G4double random = G4UniformRand();
643
644 if (particle == G4Electron::ElectronDefinition())
645 {
646 if((eVecEZ[z].at(0)<=energy)&&(energy<eVecEZ[z].back()))
647 {
648 auto k2
649 = std::upper_bound(eVecEZ[z].begin(),eVecEZ[z].end(), energy);
650 auto k1 = k2-1;
651
652 if ( random < eProbaShellMapZ[z][ishell][(*k1)].back()
653 && random < eProbaShellMapZ[z][ishell][(*k2)].back() )
654 {
655 auto xs12 =
656 std::upper_bound(eProbaShellMapZ[z][ishell][(*k1)].begin(),
657 eProbaShellMapZ[z][ishell][(*k1)].end(), random);
658 auto xs11 = xs12-1;
659
660 auto xs22 =
661 std::upper_bound(eProbaShellMapZ[z][ishell][(*k2)].begin(),
662 eProbaShellMapZ[z][ishell][(*k2)].end(), random);
663 auto xs21 = xs22-1;
664
665 valueE1 =*k1;
666 valueE2 =*k2;
667 valueXS21 =*xs21;
668 valueXS22 =*xs22;
669 valueXS12 =*xs12;
670 valueXS11 =*xs11;
671
672 ejeE11 = eEjectedEnergyDataZ[z][ishell][valueE1][valueXS11];
673 ejeE12 = eEjectedEnergyDataZ[z][ishell][valueE1][valueXS12];
674 ejeE21 = eEjectedEnergyDataZ[z][ishell][valueE2][valueXS21];
675 ejeE22 = eEjectedEnergyDataZ[z][ishell][valueE2][valueXS22];
676
677 secondaryEnergy = QuadInterpolator( valueXS11, valueXS12,
678 valueXS21, valueXS22,
679 ejeE11 , ejeE12 ,
680 ejeE21 , ejeE22 ,
681 valueE1, valueE2,
682 energy, random );
683 }
684 }
685 }
686 }
687
688 if(secondaryEnergy<0) secondaryEnergy=0;
689 return secondaryEnergy;
690 }
691
692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
693
694 G4ThreeVector G4DNARelativisticIonisationModel::GetEjectedElectronDirection(
695 const G4ParticleDefinition* ,
696 G4double energy,G4double secondaryenergy)
697 {
698 G4double phi = 2*CLHEP::pi*G4UniformRand();
699 G4double sintheta = std::sqrt((1.-secondaryenergy/energy)
700 / (1.+secondaryenergy/(2*CLHEP::electron_mass_c2)));
701
702 G4double dirX = sintheta*std::cos(phi);
703 G4double dirY = sintheta*std::sin(phi);
704 G4double dirZ = std::sqrt(1.-sintheta*sintheta);
705
706 G4ThreeVector vec(dirX,dirY,dirZ);
707 return vec;
708 }
709
710 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
711
712 G4double G4DNARelativisticIonisationModel::Interpolate( G4double e1,
713 G4double e2,
714 G4double e,
715 G4double xs1,
716 G4double xs2)
717 {
718
719 G4double value = 0.;
720
721 if((xs1!=0)&&(e1!=0)){
722 // Log-log interpolation by default
723 G4double a = (std::log10(xs2)-std::log10(xs1))
724 / (std::log10(e2)-std::log10(e1));
725 G4double b = std::log10(xs2) - a*std::log10(e2);
726 G4double sigma = a*std::log10(e) + b;
727 value = (std::pow(10.,sigma));
728 }
729 else{
730 // Lin-Lin interpolation
731 G4double d1 = xs1;
732 G4double d2 = xs2;
733 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
734 }
735
736 return value;
737 }
738
739 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
740
741 G4double G4DNARelativisticIonisationModel::QuadInterpolator(
742 G4double e11, G4double e12,
743 G4double e21, G4double e22,
744 G4double xs11, G4double xs12,
745 G4double xs21, G4double xs22,
746 G4double t1, G4double t2,
747 G4double t, G4double e)
748 {
749 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
750 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
751 G4double value
752 = Interpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
753 return value;
754 }
755
756