| 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 // Based on the work described in
27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo
29 //
30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 //
35
36 #include "G4DNAEmfietzoglouIonisationModel.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4UAtomicDeexcitation.hh"
40 #include "G4LossTableManager.hh"
41 #include "G4DNAChemistryManager.hh"
42 #include "G4DNAMolecularMaterial.hh"
43 #include "G4DNABornAngle.hh"
44 #include "G4DeltaAngle.hh"
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48 using namespace std;
49
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
52 G4DNAEmfietzoglouIonisationModel::G4DNAEmfietzoglouIonisationModel(const G4ParticleDefinition*,
53 const G4String& nam) :
54 G4VEmModel(nam)
55 {
56 verboseLevel = 0;
57 // Verbosity scale:
58 // 0 = nothing
59 // 1 = warning for energy non-conservation
60 // 2 = details of energy budget
61 // 3 = calculation of cross sections, file openings, sampling of atoms
62 // 4 = entering in methods
63
64 if(verboseLevel > 0)
65 {
66 G4cout << "Emfietzoglou ionisation model is constructed " << G4endl;
67 }
68
69 // Mark this model as "applicable" for atomic deexcitation
70 SetDeexcitationFlag(true);
71 fAtomDeexcitation = nullptr;
72 fParticleChangeForGamma = nullptr;
73 fpMolWaterDensity = nullptr;
74
75 // Define default angular generator
76 SetAngularDistribution(new G4DNABornAngle());
77
78 SetLowEnergyLimit(10. * eV);
79 SetHighEnergyLimit(10. * keV);
80
81 // Selection of computation method
82
83 fasterCode = false;
84
85 // Selection of stationary mode
86
87 statCode = false;
88 }
89
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91
92 G4DNAEmfietzoglouIonisationModel::~G4DNAEmfietzoglouIonisationModel()
93 {
94 // Cross section
95
96 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
97 for(pos = tableData.begin(); pos != tableData.end(); ++pos)
98 {
99 G4DNACrossSectionDataSet* table = pos->second;
100 delete table;
101 }
102
103 // Final state
104
105 eVecm.clear();
106
107 }
108
109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
111 void G4DNAEmfietzoglouIonisationModel::Initialise(const G4ParticleDefinition* particle,
112 const G4DataVector& /*cuts*/)
113 {
114
115 if(verboseLevel > 3)
116 {
117 G4cout << "Calling G4DNAEmfietzoglouIonisationModel::Initialise()" << G4endl;
118 }
119
120 // Energy limits
121
122 G4String fileElectron("dna/sigma_ionisation_e_emfietzoglou");
123
124 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
125
126 G4String electron;
127
128 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
129
130 const char *path = G4FindDataDir("G4LEDATA");
131
132 // *** ELECTRON
133
134 electron = electronDef->GetParticleName();
135
136 tableFile[electron] = fileElectron;
137
138 // Cross section
139
140 auto tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141 tableE->LoadData(fileElectron);
142
143 tableData[electron] = tableE;
144
145 // Final state
146
147 std::ostringstream eFullFileName;
148
149 if (fasterCode) eFullFileName << path << "/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat";
150 if (!fasterCode) eFullFileName << path << "/dna/sigmadiff_ionisation_e_emfietzoglou.dat";
151
152 std::ifstream eDiffCrossSection(eFullFileName.str().c_str());
153
154 if (!eDiffCrossSection)
155 {
156 if (fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
157 FatalException,"Missing data file:/dna/sigmadiff_cumulated_ionisation_e_emfietzoglou.dat");
158
159 if (!fasterCode) G4Exception("G4DNAEmfietzoglouIonisationModel::Initialise","em0003",
160 FatalException,"Missing data file:/dna/sigmadiff_ionisation_e_emfietzoglou.dat");
161 }
162
163 //
164
165 // Clear the arrays for re-initialization case (MT mode)
166 // March 25th, 2014 - Vaclav Stepan, Sebastien Incerti
167
168 eTdummyVec.clear();
169
170 eVecm.clear();
171
172 eProbaShellMap->clear();
173
174 eDiffCrossSectionData->clear();
175
176 eNrjTransfData->clear();
177
178 //
179
180 eTdummyVec.push_back(0.);
181 while(!eDiffCrossSection.eof())
182 {
183 G4double tDummy;
184 G4double eDummy;
185 eDiffCrossSection>>tDummy>>eDummy;
186 if (tDummy != eTdummyVec.back()) eTdummyVec.push_back(tDummy);
187 for (G4int j=0; j<5; j++)
188 {
189 eDiffCrossSection>>eDiffCrossSectionData[j][tDummy][eDummy];
190
191 if (fasterCode)
192 {
193 eNrjTransfData[j][tDummy][eDiffCrossSectionData[j][tDummy][eDummy]]=eDummy;
194 eProbaShellMap[j][tDummy].push_back(eDiffCrossSectionData[j][tDummy][eDummy]);
195 }
196
197 // SI - only if eof is not reached
198 if (!eDiffCrossSection.eof() && !fasterCode)
199 {
200 eDiffCrossSectionData[j][tDummy][eDummy]*=scaleFactor;
201 }
202
203 if (!fasterCode) eVecm[tDummy].push_back(eDummy);
204
205 }
206 }
207
208 //
209
210 if( verboseLevel>0 )
211 {
212 G4cout << "Emfietzoglou ionisation model is initialized " << G4endl
213 << "Energy range: "
214 << LowEnergyLimit() / eV << " eV - "
215 << HighEnergyLimit() / keV << " keV for "
216 << particle->GetParticleName()
217 << G4endl;
218 }
219
220 // Initialize water density pointer
221
222 fpMolWaterDensity =
223 G4DNAMolecularMaterial::Instance()->
224 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
225
226 // AD
227
228 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
229
230 if (isInitialised)
231 { return;}
232 fParticleChangeForGamma = GetParticleChangeForGamma();
233 isInitialised = true;
234 }
235
236 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237
238 G4double G4DNAEmfietzoglouIonisationModel::
239 CrossSectionPerVolume(const G4Material* material,
240 const G4ParticleDefinition* particleDefinition,
241 G4double ekin,
242 G4double,
243 G4double)
244 {
245 if(verboseLevel > 3)
246 {
247 G4cout
248 << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouIonisationModel"
249 << G4endl;
250 }
251
252 if (particleDefinition != G4Electron::ElectronDefinition()) return 0; // necessary ??
253
254 // Calculate total cross section for model
255
256 G4double sigma=0;
257
258 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
259
260 const G4String& particleName = particleDefinition->GetParticleName();
261
262 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
263 {
264 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
265 pos = tableData.find(particleName);
266
267 if (pos != tableData.end())
268 {
269 G4DNACrossSectionDataSet* table = pos->second;
270 if (table != nullptr)
271 {
272 sigma = table->FindValue(ekin);
273 }
274 }
275 else
276 {
277 G4Exception("G4DNAEmfietzoglouIonisationModel::CrossSectionPerVolume","em0002",
278 FatalException,"Model not applicable to particle type.");
279 }
280 }
281
282 if (verboseLevel > 2)
283 {
284 G4cout << "__________________________________" << G4endl;
285 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO START" << G4endl;
286 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
287 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
288 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
289 G4cout << "G4DNAEmfietzoglouIonisationModel - XS INFO END" << G4endl;
290 }
291
292 return sigma*waterDensity;
293 }
294
295 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
296
297 void G4DNAEmfietzoglouIonisationModel::
298 SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
299 const G4MaterialCutsCouple* couple,
300 const G4DynamicParticle* particle,
301 G4double,
302 G4double)
303 {
304
305 if(verboseLevel > 3)
306 {
307 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouIonisationModel"
308 << G4endl;
309 }
310
311 G4double k = particle->GetKineticEnergy();
312
313 const G4String& particleName = particle->GetDefinition()->GetParticleName();
314
315 if (k >= LowEnergyLimit() && k <= HighEnergyLimit())
316 {
317 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
318 G4double particleMass = particle->GetDefinition()->GetPDGMass();
319 G4double totalEnergy = k + particleMass;
320 G4double pSquare = k * (totalEnergy + particleMass);
321 G4double totalMomentum = std::sqrt(pSquare);
322
323 G4int ionizationShell = 0;
324
325 ionizationShell = RandomSelect(k,particleName);
326
327 G4double bindingEnergy = 0;
328 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
329
330 // SI : additional protection if tcs interpolation method is modified
331 if (k<bindingEnergy) return;
332 //
333
334 G4double secondaryKinetic=-1000*eV;
335
336 if (!fasterCode) secondaryKinetic = RandomizeEjectedElectronEnergy(particle->GetDefinition(),k,ionizationShell);
337
338 if (fasterCode)
339 secondaryKinetic = RandomizeEjectedElectronEnergyFromCumulatedDcs(particle->GetDefinition(),k,ionizationShell);
340
341 // SI - For atom. deexc. tagging - 23/05/2017
342
343 G4int Z = 8;
344
345 G4ThreeVector deltaDirection =
346 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
347 Z, ionizationShell,
348 couple->GetMaterial());
349
350 if (secondaryKinetic>0)
351 {
352 auto dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
353 fvect->push_back(dp);
354 }
355
356 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
357
358 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
359 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
360 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
361 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz);
362 finalPx /= finalMomentum;
363 finalPy /= finalMomentum;
364 finalPz /= finalMomentum;
365
366 G4ThreeVector direction;
367 direction.set(finalPx,finalPy,finalPz);
368
369 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit());
370
371 // AM: sample deexcitation
372 // here we assume that H_{2}O electronic levels are the same as Oxygen.
373 // this can be considered true with a rough 10% error in energy on K-shell,
374
375 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
376 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
377
378 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
379
380 // SI: only atomic deexcitation from K shell is considered
381 if((fAtomDeexcitation != nullptr) && ionizationShell == 4)
382 {
383 const G4AtomicShell* shell
384 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
385 secNumberInit = fvect->size();
386 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
387 secNumberFinal = fvect->size();
388
389 if(secNumberFinal > secNumberInit) {
390 for (size_t i=secNumberInit; i<secNumberFinal; ++i) {
391 //Check if there is enough residual energy
392 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
393 {
394 //Ok, this is a valid secondary: keep it
395 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
396 }
397 else
398 {
399 //Invalid secondary: not enough energy to create it!
400 //Keep its energy in the local deposit
401 delete (*fvect)[i];
402 (*fvect)[i]=nullptr;
403 }
404 }
405 }
406
407 }
408
409 //This should never happen
410 if(bindingEnergy < 0.0)
411 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
412 "em2050",FatalException,"Negative local energy deposit");
413
414 //bindingEnergy has been decreased
415 //by the amount of energy taken away by deexc. products
416 if (!statCode)
417 {
418 fParticleChangeForGamma->SetProposedKineticEnergy(scatteredEnergy);
419 fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
420 }
421 else
422 {
423 fParticleChangeForGamma->SetProposedKineticEnergy(k);
424 fParticleChangeForGamma->ProposeLocalEnergyDeposit(k-scatteredEnergy);
425 }
426 // TEST //////////////////////////
427 // if (secondaryKinetic<0) abort();
428 // if (scatteredEnergy<0) abort();
429 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
430 // if (k-scatteredEnergy<0) abort();
431 /////////////////////////////////
432
433 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
434 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eIonizedMolecule,
435 ionizationShell,
436 theIncomingTrack);
437 }
438
439 }
440
441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
442
443 G4double
444 G4DNAEmfietzoglouIonisationModel::
445 RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
446 G4double k,
447 G4int shell)
448 {
449 // G4cout << "*** SLOW computation for "
450 // << " " << particleDefinition->GetParticleName() << G4endl;
451
452 if(particleDefinition == G4Electron::ElectronDefinition())
453 {
454 G4double maximumEnergyTransfer = 0.;
455 if((k + waterStructure.IonisationEnergy(shell)) / 2. > k)
456 maximumEnergyTransfer = k;
457 else
458 maximumEnergyTransfer = (k + waterStructure.IonisationEnergy(shell))/ 2.;
459
460 // SI : original method
461 /*
462 G4double crossSectionMaximum = 0.;
463 for(G4double value=waterStructure.IonisationEnergy(shell); value<=maximumEnergyTransfer; value+=0.1*eV)
464 {
465 G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k/eV, value/eV, shell);
466 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
467 }
468 */
469
470 // SI : alternative method
471 G4double crossSectionMaximum = 0.;
472
473 G4double minEnergy = waterStructure.IonisationEnergy(shell);
474 G4double maxEnergy = maximumEnergyTransfer;
475 G4int nEnergySteps = 50;
476
477 G4double value(minEnergy);
478 G4double stpEnergy(std::pow(maxEnergy / value,
479 1. / static_cast<G4double>(nEnergySteps - 1)));
480 G4int step(nEnergySteps);
481 while(step > 0)
482 {
483 step--;
484 G4double differentialCrossSection =
485 DifferentialCrossSection(particleDefinition,
486 k / eV,
487 value / eV,
488 shell);
489 if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum =
490 differentialCrossSection;
491 value *= stpEnergy;
492 }
493 //
494
495 G4double secondaryElectronKineticEnergy = 0.;
496 do
497 {
498 secondaryElectronKineticEnergy = G4UniformRand()* (maximumEnergyTransfer-waterStructure.IonisationEnergy(shell));
499 }while(G4UniformRand()*crossSectionMaximum >
500 DifferentialCrossSection(particleDefinition, k/eV,
501 (secondaryElectronKineticEnergy+waterStructure.IonisationEnergy(shell))/eV,shell));
502
503 return secondaryElectronKineticEnergy;
504
505 }
506
507 return 0;
508 }
509
510 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
511
512 // The following section is not used anymore but is kept for memory
513 // GetAngularDistribution()->SampleDirectionForShell is used instead
514
515 /*
516 void G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
517 G4double k,
518 G4double secKinetic,
519 G4double & cosTheta,
520 G4double & phi )
521 {
522 if (particleDefinition == G4Electron::ElectronDefinition())
523 {
524 phi = twopi * G4UniformRand();
525 if (secKinetic < 50.*eV) cosTheta = (2.*G4UniformRand())-1.;
526 else if (secKinetic <= 200.*eV)
527 {
528 if (G4UniformRand() <= 0.1) cosTheta = (2.*G4UniformRand())-1.;
529 else cosTheta = G4UniformRand()*(std::sqrt(2.)/2);
530 }
531 else
532 {
533 G4double sin2O = (1.-secKinetic/k) / (1.+secKinetic/(2.*electron_mass_c2));
534 cosTheta = std::sqrt(1.-sin2O);
535 }
536 }
537
538 else if (particleDefinition == G4Proton::ProtonDefinition())
539 {
540 G4double maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
541 phi = twopi * G4UniformRand();
542
543 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
544
545 // Restriction below 100 eV from Emfietzoglou (2000)
546
547 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
548 else cosTheta = (2.*G4UniformRand())-1.;
549
550 }
551 }
552 */
553
554 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
555 G4double G4DNAEmfietzoglouIonisationModel::DifferentialCrossSection(G4ParticleDefinition * particleDefinition,
556 G4double k,
557 G4double energyTransfer,
558 G4int ionizationLevelIndex)
559 {
560 G4double sigma = 0.;
561
562 if(energyTransfer >= waterStructure.IonisationEnergy(ionizationLevelIndex)/eV)
563 {
564 G4double valueT1 = 0;
565 G4double valueT2 = 0;
566 G4double valueE21 = 0;
567 G4double valueE22 = 0;
568 G4double valueE12 = 0;
569 G4double valueE11 = 0;
570
571 G4double xs11 = 0;
572 G4double xs12 = 0;
573 G4double xs21 = 0;
574 G4double xs22 = 0;
575
576 if(particleDefinition == G4Electron::ElectronDefinition())
577 {
578 // Protection against out of boundary access
579 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
580 //
581
582 // k should be in eV and energy transfer eV also
583
584 auto t2 = std::upper_bound(eTdummyVec.begin(),
585 eTdummyVec.end(),
586 k);
587
588 auto t1 = t2 - 1;
589
590 // SI : the following condition avoids situations where energyTransfer >last vector element
591 // added strict limitations (09/08/2017)
592 if(energyTransfer < eVecm[(*t1)].back() &&
593 energyTransfer < eVecm[(*t2)].back())
594 {
595 auto e12 =
596 std::upper_bound(eVecm[(*t1)].begin(),
597 eVecm[(*t1)].end(),
598 energyTransfer);
599 auto e11 = e12 - 1;
600
601 auto e22 =
602 std::upper_bound(eVecm[(*t2)].begin(),
603 eVecm[(*t2)].end(),
604 energyTransfer);
605 auto e21 = e22 - 1;
606
607 valueT1 = *t1;
608 valueT2 = *t2;
609 valueE21 = *e21;
610 valueE22 = *e22;
611 valueE12 = *e12;
612 valueE11 = *e11;
613
614 xs11 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE11];
615 xs12 = eDiffCrossSectionData[ionizationLevelIndex][valueT1][valueE12];
616 xs21 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE21];
617 xs22 = eDiffCrossSectionData[ionizationLevelIndex][valueT2][valueE22];
618
619 //G4cout << "-------------------" << G4endl;
620 //G4cout << "ionizationLevelIndex=" << ionizationLevelIndex << G4endl;
621 //G4cout << "valueT1/eV=" << valueT1 << " valueT2/eV=" << valueT2 << G4endl;
622 //G4cout << "valueE11/eV=" << valueE11 << " valueE12/eV=" << valueE12
623 // << " valueE21/eV=" << valueE21 << " valueE22/eV=" << valueE22 << G4endl;
624 //G4cout << "xs11=" << xs11 / ((1.e-22 / 3.343) * m*m) << G4endl;
625 //G4cout << "xs12=" << xs12 / ((1.e-22 / 3.343) * m*m) << G4endl;
626 //G4cout << "xs21=" << xs21 / ((1.e-22 / 3.343) * m*m) << G4endl;
627 //G4cout << "xs22=" << xs22 / ((1.e-22 / 3.343) * m*m) << G4endl;
628 //G4cout << "###################" << G4endl;
629
630 }
631
632 }
633
634 G4double xsProduct = xs11 * xs12 * xs21 * xs22;
635 if(xsProduct != 0.)
636 {
637 sigma = QuadInterpolator(valueE11,
638 valueE12,
639 valueE21,
640 valueE22,
641 xs11,
642 xs12,
643 xs21,
644 xs22,
645 valueT1,
646 valueT2,
647 k,
648 energyTransfer);
649 }
650
651 }
652
653 return sigma;
654 }
655
656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
657
658 G4double G4DNAEmfietzoglouIonisationModel::Interpolate(G4double e1,
659 G4double e2,
660 G4double e,
661 G4double xs1,
662 G4double xs2)
663 {
664
665 G4double value = 0.;
666
667 // Log-log interpolation by default
668
669 if(e1 != 0 && e2 != 0 && (std::log10(e2) - std::log10(e1)) != 0
670 && !fasterCode)
671 {
672 G4double a = (std::log10(xs2) - std::log10(xs1))
673 / (std::log10(e2) - std::log10(e1));
674 G4double b = std::log10(xs2) - a * std::log10(e2);
675 G4double sigma = a * std::log10(e) + b;
676 value = (std::pow(10., sigma));
677 }
678
679 // Switch to lin-lin interpolation
680 /*
681 if ((e2-e1)!=0)
682 {
683 G4double d1 = xs1;
684 G4double d2 = xs2;
685 value = (d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
686 }
687 */
688
689 // Switch to log-lin interpolation for faster code
690 if((e2 - e1) != 0 && xs1 != 0 && xs2 != 0 && fasterCode)
691 {
692 G4double d1 = std::log10(xs1);
693 G4double d2 = std::log10(xs2);
694 value = std::pow(10., (d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
695 }
696
697 // Switch to lin-lin interpolation for faster code
698 // in case one of xs1 or xs2 (=cum proba) value is zero
699
700 if((e2 - e1) != 0 && (xs1 == 0 || xs2 == 0) && fasterCode)
701 {
702 G4double d1 = xs1;
703 G4double d2 = xs2;
704 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
705 }
706
707 /*
708 G4cout
709 << e1 << " "
710 << e2 << " "
711 << e << " "
712 << xs1 << " "
713 << xs2 << " "
714 << value
715 << G4endl;
716 */
717
718 return value;
719 }
720
721 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
722
723 G4double G4DNAEmfietzoglouIonisationModel::QuadInterpolator(G4double e11,
724 G4double e12,
725 G4double e21,
726 G4double e22,
727 G4double xs11,
728 G4double xs12,
729 G4double xs21,
730 G4double xs22,
731 G4double t1,
732 G4double t2,
733 G4double t,
734 G4double e)
735 {
736 G4double interpolatedvalue1 = Interpolate(e11, e12, e, xs11, xs12);
737 G4double interpolatedvalue2 = Interpolate(e21, e22, e, xs21, xs22);
738 G4double value = Interpolate(t1,
739 t2,
740 t,
741 interpolatedvalue1,
742 interpolatedvalue2);
743
744 return value;
745 }
746
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
748
749 G4int G4DNAEmfietzoglouIonisationModel::RandomSelect(G4double k,
750 const G4String& particle)
751 {
752 G4int level = 0;
753
754 auto pos = tableData.find(particle);
755
756 if(pos != tableData.cend())
757 {
758 G4DNACrossSectionDataSet* table = pos->second;
759
760 if(table != nullptr)
761 {
762 auto valuesBuffer = new G4double[table->NumberOfComponents()];
763 const auto n = (G4int)table->NumberOfComponents();
764 G4int i(n);
765 G4double value = 0.;
766
767 while(i > 0)
768 {
769 i--;
770 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
771 value += valuesBuffer[i];
772 }
773
774 value *= G4UniformRand();
775
776 i = n;
777
778 while(i > 0)
779 {
780 i--;
781
782 if(valuesBuffer[i] > value)
783 {
784 delete[] valuesBuffer;
785 return i;
786 }
787 value -= valuesBuffer[i];
788 }
789
790 delete[] valuesBuffer;
791
792 }
793 }
794 else
795 {
796 G4Exception("G4DNAEmfietzoglouIonisationModel::RandomSelect",
797 "em0002",
798 FatalException,
799 "Model not applicable to particle type.");
800 }
801
802 return level;
803 }
804
805 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
806
807 G4double G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
808 G4double k,
809 G4int shell)
810 {
811 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
812
813 G4double secondaryElectronKineticEnergy = 0.;
814
815 secondaryElectronKineticEnergy = RandomTransferedEnergy(particleDefinition,
816 k / eV,
817 shell)
818 * eV
819 - waterStructure.IonisationEnergy(shell);
820
821 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
822 if(secondaryElectronKineticEnergy < 0.) return 0.;
823
824 return secondaryElectronKineticEnergy;
825 }
826
827 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
828
829 G4double G4DNAEmfietzoglouIonisationModel::RandomTransferedEnergy(G4ParticleDefinition* particleDefinition,
830 G4double k,
831 G4int ionizationLevelIndex)
832 {
833
834 G4double random = G4UniformRand();
835
836 G4double nrj = 0.;
837
838 G4double valueK1 = 0;
839 G4double valueK2 = 0;
840 G4double valuePROB21 = 0;
841 G4double valuePROB22 = 0;
842 G4double valuePROB12 = 0;
843 G4double valuePROB11 = 0;
844
845 G4double nrjTransf11 = 0;
846 G4double nrjTransf12 = 0;
847 G4double nrjTransf21 = 0;
848 G4double nrjTransf22 = 0;
849
850 if (particleDefinition == G4Electron::ElectronDefinition())
851 {
852 // Protection against out of boundary access
853 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
854 //
855
856 // k should be in eV
857 auto k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
858
859 auto k1 = k2-1;
860
861 /*
862 G4cout << "----> k=" << k
863 << " " << *k1
864 << " " << *k2
865 << " " << random
866 << " " << ionizationLevelIndex
867 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
868 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
869 << G4endl;
870 */
871
872 // SI : the following condition avoids situations where random >last vector element
873 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
874 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
875
876 {
877 auto prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
878 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
879
880 auto prob11 = prob12-1;
881
882 auto prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
883 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
884
885 auto prob21 = prob22-1;
886
887 valueK1 =*k1;
888 valueK2 =*k2;
889 valuePROB21 =*prob21;
890 valuePROB22 =*prob22;
891 valuePROB12 =*prob12;
892 valuePROB11 =*prob11;
893
894 /*
895 G4cout << " " << random << " " << valuePROB11 << " "
896 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
897 */
898
899 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
900 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
901 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
902 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
903
904 /*
905 G4cout << " " << ionizationLevelIndex << " "
906 << random << " " <<valueK1 << " " << valueK2 << G4endl;
907
908 G4cout << " " << random << " " << nrjTransf11 << " "
909 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
910 */
911
912 }
913
914 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
915
916 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
917
918 {
919 auto prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
920 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
921
922 auto prob21 = prob22-1;
923
924 valueK1 =*k1;
925 valueK2 =*k2;
926 valuePROB21 =*prob21;
927 valuePROB22 =*prob22;
928
929 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
930
931 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
932 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
933
934 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
935
936 // zeros are explicitly set
937
938 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
939
940 /*
941 G4cout << " " << ionizationLevelIndex << " "
942 << random << " " <<valueK1 << " " << valueK2 << G4endl;
943
944 G4cout << " " << random << " " << nrjTransf11 << " "
945 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
946
947 G4cout << "ici" << " " << value << G4endl;
948 */
949
950 return value;
951 }
952
953 }
954
955 // End electron
956
957 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
958
959 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
960
961 if (nrjTransfProduct != 0.)
962 {
963 nrj = QuadInterpolator( valuePROB11, valuePROB12,
964 valuePROB21, valuePROB22,
965 nrjTransf11, nrjTransf12,
966 nrjTransf21, nrjTransf22,
967 valueK1, valueK2,
968 k, random);
969 }
970
971 //G4cout << nrj << endl;
972
973 return nrj;
974 }
975