Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4DNAEmfietzoglouIonisationModel.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// 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 ([email protected])
32// D. Emfietzoglou ([email protected])
33// S. Incerti ([email protected])
34//
35
38#include "G4SystemOfUnits.hh"
40#include "G4LossTableManager.hh"
43#include "G4DNABornAngle.hh"
44#include "G4DeltaAngle.hh"
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48using namespace std;
49
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
51
53 const G4String& nam) :
54G4VEmModel(nam), isInitialised(false)
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
71 fAtomDeexcitation = 0;
73 fpMolWaterDensity = 0;
74
75 // Define default angular generator
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
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
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
125
126 G4String electron;
127
128 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
129
130 char *path = getenv("G4LEDATA");
131
132 // *** ELECTRON
133
134 electron = electronDef->GetParticleName();
135
136 tableFile[electron] = fileElectron;
137
138 // Cross section
139
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 =
224 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
225
226 // AD
227
228 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
229
230 if (isInitialised)
231 { return;}
233 isInitialised = true;
234}
235
236//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
237
239CrossSectionPerVolume(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 != 0)
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
298SampleSecondaries(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 G4DynamicParticle* 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
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 && 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]=0;
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 {
420 }
421 else
422 {
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();
435 ionizationShell,
436 theIncomingTrack);
437 }
438
439}
440
441//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
442
444G4DNAEmfietzoglouIonisationModel::
445RandomizeEjectedElectronEnergy(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......
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 std::vector<G4double>::iterator t2 = std::upper_bound(eTdummyVec.begin(),
585 eTdummyVec.end(),
586 k);
587
588 std::vector<G4double>::iterator 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 std::vector<G4double>::iterator e12 =
596 std::upper_bound(eVecm[(*t1)].begin(),
597 eVecm[(*t1)].end(),
598 energyTransfer);
599 std::vector<G4double>::iterator e11 = e12 - 1;
600
601 std::vector<G4double>::iterator e22 =
602 std::upper_bound(eVecm[(*t2)].begin(),
603 eVecm[(*t2)].end(),
604 energyTransfer);
605 std::vector<G4double>::iterator 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
658G4double 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
723G4double 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
749G4int G4DNAEmfietzoglouIonisationModel::RandomSelect(G4double k,
750 const G4String& particle)
751{
752 G4int level = 0;
753
754 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
755 pos = tableData.find(particle);
756
757 if(pos != tableData.end())
758 {
759 G4DNACrossSectionDataSet* table = pos->second;
760
761 if(table != 0)
762 {
763 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
764 const size_t n(table->NumberOfComponents());
765 size_t i(n);
766 G4double value = 0.;
767
768 while(i > 0)
769 {
770 i--;
771 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
772 value += valuesBuffer[i];
773 }
774
775 value *= G4UniformRand();
776
777 i = n;
778
779 while(i > 0)
780 {
781 i--;
782
783 if(valuesBuffer[i] > value)
784 {
785 delete[] valuesBuffer;
786 return i;
787 }
788 value -= valuesBuffer[i];
789 }
790
791 if(valuesBuffer) delete[] valuesBuffer;
792
793 }
794 }
795 else
796 {
797 G4Exception("G4DNAEmfietzoglouIonisationModel::RandomSelect",
798 "em0002",
800 "Model not applicable to particle type.");
801 }
802
803 return level;
804}
805
806//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
807
808G4double G4DNAEmfietzoglouIonisationModel::RandomizeEjectedElectronEnergyFromCumulatedDcs(G4ParticleDefinition* particleDefinition,
809 G4double k,
810 G4int shell)
811{
812 //G4cout << "*** FAST computation for " << " " << particleDefinition->GetParticleName() << G4endl;
813
814 G4double secondaryElectronKineticEnergy = 0.;
815
816 secondaryElectronKineticEnergy = RandomTransferedEnergy(particleDefinition,
817 k / eV,
818 shell)
819 * eV
820 - waterStructure.IonisationEnergy(shell);
821
822 //G4cout << RandomTransferedEnergy(particleDefinition, k/eV, shell) << G4endl;
823 if(secondaryElectronKineticEnergy < 0.) return 0.;
824
825 return secondaryElectronKineticEnergy;
826}
827
828//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
829
830G4double G4DNAEmfietzoglouIonisationModel::RandomTransferedEnergy(G4ParticleDefinition* particleDefinition,
831 G4double k,
832 G4int ionizationLevelIndex)
833{
834
835 G4double random = G4UniformRand();
836
837 G4double nrj = 0.;
838
839 G4double valueK1 = 0;
840 G4double valueK2 = 0;
841 G4double valuePROB21 = 0;
842 G4double valuePROB22 = 0;
843 G4double valuePROB12 = 0;
844 G4double valuePROB11 = 0;
845
846 G4double nrjTransf11 = 0;
847 G4double nrjTransf12 = 0;
848 G4double nrjTransf21 = 0;
849 G4double nrjTransf22 = 0;
850
851 if (particleDefinition == G4Electron::ElectronDefinition())
852 {
853 // Protection against out of boundary access
854 if (k==eTdummyVec.back()) k=k*(1.-1e-12);
855 //
856
857 // k should be in eV
858 std::vector<G4double>::iterator k2 = std::upper_bound(eTdummyVec.begin(),eTdummyVec.end(), k);
859
860 std::vector<G4double>::iterator k1 = k2-1;
861
862 /*
863 G4cout << "----> k=" << k
864 << " " << *k1
865 << " " << *k2
866 << " " << random
867 << " " << ionizationLevelIndex
868 << " " << eProbaShellMap[ionizationLevelIndex][(*k1)].back()
869 << " " << eProbaShellMap[ionizationLevelIndex][(*k2)].back()
870 << G4endl;
871 */
872
873 // SI : the following condition avoids situations where random >last vector element
874 if ( random <= eProbaShellMap[ionizationLevelIndex][(*k1)].back()
875 && random <= eProbaShellMap[ionizationLevelIndex][(*k2)].back() )
876
877 {
878 std::vector<G4double>::iterator prob12 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k1)].begin(),
879 eProbaShellMap[ionizationLevelIndex][(*k1)].end(), random);
880
881 std::vector<G4double>::iterator prob11 = prob12-1;
882
883 std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
884 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
885
886 std::vector<G4double>::iterator prob21 = prob22-1;
887
888 valueK1 =*k1;
889 valueK2 =*k2;
890 valuePROB21 =*prob21;
891 valuePROB22 =*prob22;
892 valuePROB12 =*prob12;
893 valuePROB11 =*prob11;
894
895 /*
896 G4cout << " " << random << " " << valuePROB11 << " "
897 << valuePROB12 << " " << valuePROB21 << " " << valuePROB22 << G4endl;
898 */
899
900 nrjTransf11 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB11];
901 nrjTransf12 = eNrjTransfData[ionizationLevelIndex][valueK1][valuePROB12];
902 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
903 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
904
905 /*
906 G4cout << " " << ionizationLevelIndex << " "
907 << random << " " <<valueK1 << " " << valueK2 << G4endl;
908
909 G4cout << " " << random << " " << nrjTransf11 << " "
910 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
911 */
912
913 }
914
915 // Avoids cases where cum xs is zero for k1 and is not for k2 (with always k1<k2)
916
917 if ( random > eProbaShellMap[ionizationLevelIndex][(*k1)].back() )
918
919 {
920 std::vector<G4double>::iterator prob22 = std::upper_bound(eProbaShellMap[ionizationLevelIndex][(*k2)].begin(),
921 eProbaShellMap[ionizationLevelIndex][(*k2)].end(), random);
922
923 std::vector<G4double>::iterator prob21 = prob22-1;
924
925 valueK1 =*k1;
926 valueK2 =*k2;
927 valuePROB21 =*prob21;
928 valuePROB22 =*prob22;
929
930 //G4cout << " " << random << " " << valuePROB21 << " " << valuePROB22 << G4endl;
931
932 nrjTransf21 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB21];
933 nrjTransf22 = eNrjTransfData[ionizationLevelIndex][valueK2][valuePROB22];
934
935 G4double interpolatedvalue2 = Interpolate(valuePROB21, valuePROB22, random, nrjTransf21, nrjTransf22);
936
937 // zeros are explicitly set
938
939 G4double value = Interpolate(valueK1, valueK2, k, 0., interpolatedvalue2);
940
941 /*
942 G4cout << " " << ionizationLevelIndex << " "
943 << random << " " <<valueK1 << " " << valueK2 << G4endl;
944
945 G4cout << " " << random << " " << nrjTransf11 << " "
946 << nrjTransf12 << " " << nrjTransf21 << " " <<nrjTransf22 << G4endl;
947
948 G4cout << "ici" << " " << value << G4endl;
949 */
950
951 return value;
952 }
953
954 }
955
956 // End electron
957
958 G4double nrjTransfProduct = nrjTransf11 * nrjTransf12 * nrjTransf21 * nrjTransf22;
959
960 //G4cout << "nrjTransfProduct=" << nrjTransfProduct << G4endl;
961
962 if (nrjTransfProduct != 0.)
963 {
964 nrj = QuadInterpolator( valuePROB11, valuePROB12,
965 valuePROB21, valuePROB22,
966 nrjTransf11, nrjTransf12,
967 nrjTransf21, nrjTransf22,
968 valueK1, valueK2,
969 k, random);
970 }
971
972 //G4cout << nrj << endl;
973
974 return nrj;
975}
G4AtomicShellEnumerator
@ eIonizedMolecule
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
double z() const
Hep3Vector unit() const
double x() const
double y() const
void set(double x, double y, double z)
static G4DNAChemistryManager * Instance()
void CreateWaterMolecule(ElectronicModification, G4int, const G4Track *)
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual size_t NumberOfComponents(void) const
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual G4bool LoadData(const G4String &argFileName)
G4DNAEmfietzoglouIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNAEmfietzoglouIonisationModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &= *(new G4DataVector()))
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
G4double DifferentialCrossSection(G4ParticleDefinition *aParticleDefinition, G4double k, G4double energyTransfer, G4int shell)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static G4DNAMolecularMaterial * Instance()
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:88
static G4Electron * Electron()
Definition: G4Electron.cc:93
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
size_t GetIndex() const
Definition: G4Material.hh:258
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:651
const G4Track * GetCurrentTrack() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
const G4String & GetParticleName() const
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4double FindValue(G4double x, G4int componentId=0) const =0
virtual G4ThreeVector & SampleDirectionForShell(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, G4int shellID, const G4Material *)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:757
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:611
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:133
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:645
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:764
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:813
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:618
void ProposeLocalEnergyDeposit(G4double anEnergyPart)