Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4DNARuddIonisationModel.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//
27
30#include "G4SystemOfUnits.hh"
32#include "G4LossTableManager.hh"
35#include "G4DNARuddAngle.hh"
36#include "G4DeltaAngle.hh"
37#include "G4Exp.hh"
38#include "G4Pow.hh"
39#include "G4Log.hh"
40#include "G4Alpha.hh"
41
42static G4Pow * gpow = G4Pow::GetInstance();
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
45using namespace std;
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
50 const G4String& nam) :
51G4VEmModel(nam), isInitialised(false)
52{
53 slaterEffectiveCharge[0] = 0.;
54 slaterEffectiveCharge[1] = 0.;
55 slaterEffectiveCharge[2] = 0.;
56 sCoefficient[0] = 0.;
57 sCoefficient[1] = 0.;
58 sCoefficient[2] = 0.;
59
60 lowEnergyLimitForZ1 = 0 * eV;
61 lowEnergyLimitForZ2 = 0 * eV;
62 lowEnergyLimitOfModelForZ1 = 100 * eV;
63 lowEnergyLimitOfModelForZ2 = 1 * keV;
64 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
65 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
66
67 verboseLevel = 0;
68 // Verbosity scale:
69 // 0 = nothing
70 // 1 = warning for energy non-conservation
71 // 2 = details of energy budget
72 // 3 = calculation of cross sections, file openings, sampling of atoms
73 // 4 = entering in methods
74
75 if (verboseLevel > 0)
76 {
77 G4cout << "Rudd ionisation model is constructed " << G4endl;
78 }
79
80 // Define default angular generator
82
83 // Mark this model as "applicable" for atomic deexcitation
85
86 // Selection of stationary mode
87
88 statCode = false;
89}
90
91//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92
94{
95 // Cross section
96
97 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
98 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
99 {
100 G4DNACrossSectionDataSet* table = pos->second;
101 delete table;
102 }
103
104 // The following removal is forbidden since G4VEnergyLossmodel takes care of deletion
105 // Coverity however will signal this as an error
106 // if (fAtomDeexcitation) {delete fAtomDeexcitation;}
107
108}
109
110//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111
113 const G4DataVector& /*cuts*/)
114{
115
116 if (verboseLevel > 3)
117 {
118 G4cout << "Calling G4DNARuddIonisationModel::Initialise()" << G4endl;
119 }
120
121 // Energy limits
122
123 G4String fileProton("dna/sigma_ionisation_p_rudd");
124 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
125 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
126 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
127 G4String fileHelium("dna/sigma_ionisation_he_rudd");
128
129 G4DNAGenericIonsManager *instance;
131 protonDef = G4Proton::ProtonDefinition();
132 hydrogenDef = instance->GetIon("hydrogen");
133 alphaPlusPlusDef = G4Alpha::Alpha();
134 alphaPlusDef = instance->GetIon("alpha+");
135 heliumDef = instance->GetIon("helium");
136
137 G4String proton;
138 G4String hydrogen;
139 G4String alphaPlusPlus;
140 G4String alphaPlus;
141 G4String helium;
142
143 G4double scaleFactor = 1 * m*m;
144
145 // LIMITS AND DATA
146
147 // ********************************************************
148
149 proton = protonDef->GetParticleName();
150 tableFile[proton] = fileProton;
151
152 lowEnergyLimit[proton] = lowEnergyLimitForZ1;
153 highEnergyLimit[proton] = 500. * keV;
154
155 // Cross section
156
158 eV,
159 scaleFactor );
160 tableProton->LoadData(fileProton);
161 tableData[proton] = tableProton;
162
163 // ********************************************************
164
165 hydrogen = hydrogenDef->GetParticleName();
166 tableFile[hydrogen] = fileHydrogen;
167
168 lowEnergyLimit[hydrogen] = lowEnergyLimitForZ1;
169 highEnergyLimit[hydrogen] = 100. * MeV;
170
171 // Cross section
172
174 eV,
175 scaleFactor );
176 tableHydrogen->LoadData(fileHydrogen);
177
178 tableData[hydrogen] = tableHydrogen;
179
180 // ********************************************************
181
182 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
183 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
184
185 lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
186 highEnergyLimit[alphaPlusPlus] = 400. * MeV;
187
188 // Cross section
189
191 eV,
192 scaleFactor );
193 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
194
195 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
196
197 // ********************************************************
198
199 alphaPlus = alphaPlusDef->GetParticleName();
200 tableFile[alphaPlus] = fileAlphaPlus;
201
202 lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
203 highEnergyLimit[alphaPlus] = 400. * MeV;
204
205 // Cross section
206
208 eV,
209 scaleFactor );
210 tableAlphaPlus->LoadData(fileAlphaPlus);
211 tableData[alphaPlus] = tableAlphaPlus;
212
213 // ********************************************************
214
215 helium = heliumDef->GetParticleName();
216 tableFile[helium] = fileHelium;
217
218 lowEnergyLimit[helium] = lowEnergyLimitForZ2;
219 highEnergyLimit[helium] = 400. * MeV;
220
221 // Cross section
222
224 eV,
225 scaleFactor );
226 tableHelium->LoadData(fileHelium);
227 tableData[helium] = tableHelium;
228
229 //
230
231 if (particle==protonDef)
232 {
233 SetLowEnergyLimit(lowEnergyLimit[proton]);
234 SetHighEnergyLimit(highEnergyLimit[proton]);
235 }
236
237 if (particle==hydrogenDef)
238 {
239 SetLowEnergyLimit(lowEnergyLimit[hydrogen]);
240 SetHighEnergyLimit(highEnergyLimit[hydrogen]);
241 }
242
243 if (particle==heliumDef)
244 {
245 SetLowEnergyLimit(lowEnergyLimit[helium]);
246 SetHighEnergyLimit(highEnergyLimit[helium]);
247 }
248
249 if (particle==alphaPlusDef)
250 {
251 SetLowEnergyLimit(lowEnergyLimit[alphaPlus]);
252 SetHighEnergyLimit(highEnergyLimit[alphaPlus]);
253 }
254
255 if (particle==alphaPlusPlusDef)
256 {
257 SetLowEnergyLimit(lowEnergyLimit[alphaPlusPlus]);
258 SetHighEnergyLimit(highEnergyLimit[alphaPlusPlus]);
259 }
260
261 if( verboseLevel>0 )
262 {
263 G4cout << "Rudd ionisation model is initialized " << G4endl
264 << "Energy range: "
265 << LowEnergyLimit() / eV << " eV - "
266 << HighEnergyLimit() / keV << " keV for "
267 << particle->GetParticleName()
268 << G4endl;
269 }
270
271 // Initialize water density pointer
273
274 //
275
276 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
277
278 if (isInitialised)
279 { return;}
281 isInitialised = true;
282
283}
284
285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
286
288 const G4ParticleDefinition* particleDefinition,
289 G4double k,
290 G4double,
291 G4double)
292{
293 if (verboseLevel > 3)
294 {
295 G4cout << "Calling CrossSectionPerVolume() of G4DNARuddIonisationModel"
296 << G4endl;
297 }
298
299 // Calculate total cross section for model
300
301 if (
302 particleDefinition != protonDef
303 &&
304 particleDefinition != hydrogenDef
305 &&
306 particleDefinition != alphaPlusPlusDef
307 &&
308 particleDefinition != alphaPlusDef
309 &&
310 particleDefinition != heliumDef
311 )
312
313 return 0;
314
315 G4double lowLim = 0;
316
317 if ( particleDefinition == protonDef
318 || particleDefinition == hydrogenDef
319 )
320
321 lowLim = lowEnergyLimitOfModelForZ1;
322
323 if ( particleDefinition == alphaPlusPlusDef
324 || particleDefinition == alphaPlusDef
325 || particleDefinition == heliumDef
326 )
327
328 lowLim = lowEnergyLimitOfModelForZ2;
329
330 G4double highLim = 0;
331 G4double sigma=0;
332
333 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
334
335 const G4String& particleName = particleDefinition->GetParticleName();
336
337 // SI - the following is useless since lowLim is already defined
338 /*
339 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
340 pos1 = lowEnergyLimit.find(particleName);
341
342 if (pos1 != lowEnergyLimit.end())
343 {
344 lowLim = pos1->second;
345 }
346 */
347
348 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
349 pos2 = highEnergyLimit.find(particleName);
350
351 if (pos2 != highEnergyLimit.end())
352 {
353 highLim = pos2->second;
354 }
355
356 if (k <= highLim)
357 {
358 //SI : XS must not be zero otherwise sampling of secondaries method ignored
359
360 if (k < lowLim) k = lowLim;
361
362 //
363
364 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
365 pos = tableData.find(particleName);
366
367 if (pos != tableData.end())
368 {
369 G4DNACrossSectionDataSet* table = pos->second;
370 if (table != 0)
371 {
372 sigma = table->FindValue(k);
373 }
374 }
375 else
376 {
377 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
378 FatalException,"Model not applicable to particle type.");
379 }
380
381 }
382
383 if (verboseLevel > 2)
384 {
385 G4cout << "__________________________________" << G4endl;
386 G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
387 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
388 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
389 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
390 //G4cout << " - Cross section per water molecule (cm^-1)="
391 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
392 G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
393 }
394
395 return sigma*waterDensity;
396
397}
398
399//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
400
401void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
402 const G4MaterialCutsCouple* couple,
403 const G4DynamicParticle* particle,
404 G4double,
405 G4double)
406{
407 if (verboseLevel > 3)
408 {
409 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
410 << G4endl;
411 }
412
413 G4double lowLim = 0;
414 G4double highLim = 0;
415
416 if ( particle->GetDefinition() == protonDef
417 || particle->GetDefinition() == hydrogenDef
418 )
419
420 lowLim = killBelowEnergyForZ1;
421
422 if ( particle->GetDefinition() == alphaPlusPlusDef
423 || particle->GetDefinition() == alphaPlusDef
424 || particle->GetDefinition() == heliumDef
425 )
426
427 lowLim = killBelowEnergyForZ2;
428
429 G4double k = particle->GetKineticEnergy();
430
431 const G4String& particleName = particle->GetDefinition()->GetParticleName();
432
433 // SI - the following is useless since lowLim is already defined
434 /*
435 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
436 pos1 = lowEnergyLimit.find(particleName);
437
438 if (pos1 != lowEnergyLimit.end())
439 {
440 lowLim = pos1->second;
441 }
442 */
443
444 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
445 pos2 = highEnergyLimit.find(particleName);
446
447 if (pos2 != highEnergyLimit.end())
448 {
449 highLim = pos2->second;
450 }
451
452 if (k >= lowLim && k <= highLim)
453 {
454 G4ParticleDefinition* definition = particle->GetDefinition();
455 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
456 /*
457 G4double particleMass = definition->GetPDGMass();
458 G4double totalEnergy = k + particleMass;
459 G4double pSquare = k*(totalEnergy+particleMass);
460 G4double totalMomentum = std::sqrt(pSquare);
461 */
462
463 G4int ionizationShell = RandomSelect(k,particleName);
464
465 G4double bindingEnergy = 0;
466 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
467
468 //SI: additional protection if tcs interpolation method is modified
469 if (k<bindingEnergy) return;
470 //
471
472 // SI - For atom. deexc. tagging - 23/05/2017
473 G4int Z = 8;
474 //
475
476 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
477
478 G4ThreeVector deltaDirection =
479 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
480 Z, ionizationShell,
481 couple->GetMaterial());
482
483 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
484 fvect->push_back(dp);
485
486 // Ignored for ions on electrons
487 /*
488 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
489
490 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
491 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
492 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
493 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
494 finalPx /= finalMomentum;
495 finalPy /= finalMomentum;
496 finalPz /= finalMomentum;
497
498 G4ThreeVector direction;
499 direction.set(finalPx,finalPy,finalPz);
500
501 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
502 */
503
505
506 // sample deexcitation
507 // here we assume that H_{2}O electronic levels are the same of Oxigen.
508 // this can be considered true with a rough 10% error in energy on K-shell,
509
510 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
511 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
512
513 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
514
515 // SI: only atomic deexcitation from K shell is considered
516 if(fAtomDeexcitation && ionizationShell == 4)
517 {
518 const G4AtomicShell* shell
519 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
520 secNumberInit = fvect->size();
521 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
522 secNumberFinal = fvect->size();
523
524 if(secNumberFinal > secNumberInit)
525 {
526 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
527 {
528 //Check if there is enough residual energy
529 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
530 {
531 //Ok, this is a valid secondary: keep it
532 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
533 }
534 else
535 {
536 //Invalid secondary: not enough energy to create it!
537 //Keep its energy in the local deposit
538 delete (*fvect)[i];
539 (*fvect)[i]=0;
540 }
541 }
542 }
543
544 }
545
546 //This should never happen
547 if(bindingEnergy < 0.0)
548 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
549 "em2050",FatalException,"Negative local energy deposit");
550
551 //bindingEnergy has been decreased
552 //by the amount of energy taken away by deexc. products
553 if (!statCode)
554 {
557 }
558 else
559 {
562 }
563
564 // debug
565 // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
566 // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
567 // = bindingEnergy-deexSecEnergy
568 // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
569
570 // TEST //////////////////////////
571 // if (secondaryKinetic<0) abort();
572 // if (scatteredEnergy<0) abort();
573 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
574 // if (k-scatteredEnergy<0) abort();
575 /////////////////////////////////
576
577 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
579 ionizationShell,
580 theIncomingTrack);
581 }
582
583 // SI - not useful since low energy of model is 0 eV
584
585 if (k < lowLim)
586 {
590 }
591
592}
593
594//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
595
596G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
597 G4double k,
598 G4int shell)
599{
600 G4double maximumKineticEnergyTransfer = 0.;
601
602 if (particleDefinition == protonDef
603 || particleDefinition == hydrogenDef)
604 {
605 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
606 }
607
608 else if (particleDefinition == heliumDef
609 || particleDefinition == alphaPlusDef
610 || particleDefinition == alphaPlusPlusDef)
611 {
612 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
613 }
614
615 G4double crossSectionMaximum = 0.;
616
617 for (G4double value = waterStructure.IonisationEnergy(shell);
618 value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
619 value += 0.1 * eV)
620 {
621 G4double differentialCrossSection =
622 DifferentialCrossSection(particleDefinition, k, value, shell);
623 if (differentialCrossSection >= crossSectionMaximum)
624 crossSectionMaximum = differentialCrossSection;
625 }
626
627 G4double secElecKinetic = 0.;
628
629 do
630 {
631 secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
632 }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
633 k,
634 secElecKinetic+waterStructure.IonisationEnergy(shell),
635 shell));
636
637 return (secElecKinetic);
638}
639
640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641
642// The following section is not used anymore but is kept for memory
643// GetAngularDistribution()->SampleDirectionForShell is used instead
644
645/*
646 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
647 G4double k,
648 G4double secKinetic,
649 G4double & cosTheta,
650 G4double & phi )
651 {
652 G4DNAGenericIonsManager *instance;
653 instance = G4DNAGenericIonsManager::Instance();
654
655 G4double maxSecKinetic = 0.;
656
657 if (particleDefinition == protonDef
658 || particleDefinition == hydrogenDef)
659 {
660 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
661 }
662
663 else if (particleDefinition == heliumDef
664 || particleDefinition == alphaPlusDef
665 || particleDefinition == alphaPlusPlusDef)
666 {
667 maxSecKinetic = 4.* (0.511 / 3728) * k;
668 }
669
670 phi = twopi * G4UniformRand();
671
672 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
673
674 // Restriction below 100 eV from Emfietzoglou (2000)
675
676 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
677 else cosTheta = (2.*G4UniformRand())-1.;
678
679 }
680 */
681//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
682G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
683 G4double k,
684 G4double energyTransfer,
685 G4int ionizationLevelIndex)
686{
687 // Shells ids are 0 1 2 3 4 (4 is k shell)
688 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
689 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
690 //
691 // ds S F1(nu) + w * F2(nu)
692 // ---- = G(k) * ---- -------------------------------------------
693 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
694 //
695 // w is the secondary electron kinetic Energy in eV
696 //
697 // All the other parameters can be found in Rudd's Papers
698 //
699 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
700 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
701 //
702
703 const G4int j = ionizationLevelIndex;
704
705 G4double A1;
706 G4double B1;
707 G4double C1;
708 G4double D1;
709 G4double E1;
710 G4double A2;
711 G4double B2;
712 G4double C2;
713 G4double D2;
714 G4double alphaConst;
715
716 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
717 // The following values are provided by M. dingfelder (priv. comm)
718 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
719 * eV };
720
721 if (j == 4)
722 {
723 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
724 A1 = 1.25;
725 B1 = 0.5;
726 C1 = 1.00;
727 D1 = 1.00;
728 E1 = 3.00;
729 A2 = 1.10;
730 B2 = 1.30;
731 C2 = 1.00;
732 D2 = 0.00;
733 alphaConst = 0.66;
734 } else
735 {
736 //Data For Liquid Water from Dingfelder (Protons in Water)
737 A1 = 1.02;
738 B1 = 82.0;
739 C1 = 0.45;
740 D1 = -0.80;
741 E1 = 0.38;
742 A2 = 1.07;
743 // Value provided by M. Dingfelder (priv. comm)
744 B2 = 11.6;
745 //
746 C2 = 0.60;
747 D2 = 0.04;
748 alphaConst = 0.64;
749 }
750
751 const G4double n = 2.;
752 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
753
754 G4double wBig = (energyTransfer
755 - waterStructure.IonisationEnergy(ionizationLevelIndex));
756 if (wBig < 0)
757 return 0.;
758
759 G4double w = wBig / Bj[ionizationLevelIndex];
760 // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
761 if (j == 4)
762 w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
763
764 G4double Ry = 13.6 * eV;
765
766 G4double tau = 0.;
767
768 G4bool isProtonOrHydrogen = false;
769 G4bool isHelium = false;
770
771 if (particleDefinition == protonDef
772 || particleDefinition == hydrogenDef)
773 {
774 isProtonOrHydrogen = true;
775 tau = (electron_mass_c2 / proton_mass_c2) * k;
776 }
777
778 else if (particleDefinition == heliumDef
779 || particleDefinition == alphaPlusDef
780 || particleDefinition == alphaPlusPlusDef)
781 {
782 isHelium = true;
783 tau = (0.511 / 3728.) * k;
784 }
785
786 G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
787 * gpow->powN((Ry / Bj[ionizationLevelIndex]), 2);
788 if (j == 4)
789 S = 4. * pi * Bohr_radius * Bohr_radius * n
790 * gpow->powN((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
791 2);
792
793 G4double v2 = tau / Bj[ionizationLevelIndex];
794 if (j == 4)
795 v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
796
797 G4double v = std::sqrt(v2);
798 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
799 if (j == 4)
800 wc = 4. * v2 - 2. * v
801 - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
802
803 G4double L1 = (C1 * gpow->powA(v, (D1))) / (1. + E1 * gpow->powA(v, (D1 + 4.)));
804 G4double L2 = C2 * gpow->powA(v, (D2));
805 G4double H1 = (A1 * G4Log(1. + v2)) / (v2 + (B1 / v2));
806 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
807
808 G4double F1 = L1 + H1;
809 G4double F2 = (L2 * H2) / (L2 + H2);
810
811 G4double sigma =
812 CorrectionFactor(particleDefinition, k) * Gj[j]
813 * (S / Bj[ionizationLevelIndex])
814 * ((F1 + w * F2)
815 / (gpow->powN((1. + w), 3)
816 * (1. + G4Exp(alphaConst * (w - wc) / v))));
817
818 if (j == 4)
819 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
820 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
821 * ((F1 + w * F2)
822 / (gpow->powN((1. + w), 3)
823 * (1. + G4Exp(alphaConst * (w - wc) / v))));
824
825 if ((particleDefinition == hydrogenDef)
826 && (ionizationLevelIndex == 4))
827 {
828 // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
829 sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
830 * ((F1 + w * F2)
831 / (gpow->powN((1. + w), 3)
832 * (1. + G4Exp(alphaConst * (w - wc) / v))));
833 }
834
835 // if ( particleDefinition == protonDef
836 // || particleDefinition == hydrogenDef
837 // )
838
839 if (isProtonOrHydrogen)
840 {
841 return (sigma);
842 }
843
844 if (particleDefinition == alphaPlusPlusDef)
845 {
846 slaterEffectiveCharge[0] = 0.;
847 slaterEffectiveCharge[1] = 0.;
848 slaterEffectiveCharge[2] = 0.;
849 sCoefficient[0] = 0.;
850 sCoefficient[1] = 0.;
851 sCoefficient[2] = 0.;
852 }
853
854 else if (particleDefinition == alphaPlusDef)
855 {
856 slaterEffectiveCharge[0] = 2.0;
857 // The following values are provided by M. Dingfelder (priv. comm)
858 slaterEffectiveCharge[1] = 2.0;
859 slaterEffectiveCharge[2] = 2.0;
860 //
861 sCoefficient[0] = 0.7;
862 sCoefficient[1] = 0.15;
863 sCoefficient[2] = 0.15;
864 }
865
866 else if (particleDefinition == heliumDef)
867 {
868 slaterEffectiveCharge[0] = 1.7;
869 slaterEffectiveCharge[1] = 1.15;
870 slaterEffectiveCharge[2] = 1.15;
871 sCoefficient[0] = 0.5;
872 sCoefficient[1] = 0.25;
873 sCoefficient[2] = 0.25;
874 }
875
876 // if ( particleDefinition == heliumDef
877 // || particleDefinition == alphaPlusDef
878 // || particleDefinition == alphaPlusPlusDef
879 // )
880 if (isHelium)
881 {
882 sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
883 * ((F1 + w * F2)
884 / (gpow->powN((1. + w), 3)
885 * (1. + G4Exp(alphaConst * (w - wc) / v))));
886
887 if (j == 4)
888 sigma = Gj[j]
889 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
890 * ((F1 + w * F2)
891 / (gpow->powN((1. + w), 3)
892 * (1. + G4Exp(alphaConst * (w - wc) / v))));
893
894 G4double zEff = particleDefinition->GetPDGCharge() / eplus
895 + particleDefinition->GetLeptonNumber();
896
897 zEff -= (sCoefficient[0]
898 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
899 + sCoefficient[1]
900 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
901 + sCoefficient[2]
902 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
903
904 return zEff * zEff * sigma;
905 }
906
907 return 0;
908}
909
910//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
911
912G4double G4DNARuddIonisationModel::S_1s(G4double t,
913 G4double energyTransferred,
914 G4double slaterEffectiveChg,
915 G4double shellNumber)
916{
917 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
918 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
919
920 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
921 G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
922
923 return value;
924}
925
926//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
927
928G4double G4DNARuddIonisationModel::S_2s(G4double t,
929 G4double energyTransferred,
930 G4double slaterEffectiveChg,
931 G4double shellNumber)
932{
933 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
934 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
935
936 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
937 G4double value = 1.
938 - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
939
940 return value;
941
942}
943
944//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
945
946G4double G4DNARuddIonisationModel::S_2p(G4double t,
947 G4double energyTransferred,
948 G4double slaterEffectiveChg,
949 G4double shellNumber)
950{
951 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
952 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
953
954 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
955 G4double value = 1.
956 - G4Exp(-2 * r)
957 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
958
959 return value;
960}
961
962//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
963
964G4double G4DNARuddIonisationModel::R(G4double t,
965 G4double energyTransferred,
966 G4double slaterEffectiveChg,
967 G4double shellNumber)
968{
969 // tElectron = m_electron / m_alpha * t
970 // Dingfelder, in Chattanooga 2005 proceedings, p 4
971
972 G4double tElectron = 0.511 / 3728. * t;
973 // The following values are provided by M. Dingfelder (priv. comm)
974 G4double H = 2. * 13.60569172 * eV;
975 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
976 * (slaterEffectiveChg / shellNumber);
977
978 return value;
979}
980
981//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
982
983G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
984 G4double k)
985{
986
987 if (particleDefinition == G4Proton::Proton())
988 {
989 return (1.);
990 } else if (particleDefinition == hydrogenDef)
991 {
992 G4double value = (G4Log(k / eV)/gpow->logZ(10) - 4.2) / 0.5;
993 // The following values are provided by M. Dingfelder (priv. comm)
994 return ((0.6 / (1 + G4Exp(value))) + 0.9);
995 } else
996 {
997 return (1.);
998 }
999}
1000
1001//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1002
1003G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1004 const G4String& particle)
1005{
1006
1007 // BEGIN PART 1/2 OF ELECTRON CORRECTION
1008
1009 // add ONE or TWO electron-water ionisation for alpha+ and helium
1010
1011 G4int level = 0;
1012
1013 // Retrieve data table corresponding to the current particle type
1014
1015 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1016 pos = tableData.find(particle);
1017
1018 if (pos != tableData.end())
1019 {
1020 G4DNACrossSectionDataSet* table = pos->second;
1021
1022 if (table != 0)
1023 {
1024 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1025
1026 const G4int n = (G4int)table->NumberOfComponents();
1027 G4int i(n);
1028 G4double value = 0.;
1029
1030 while (i > 0)
1031 {
1032 --i;
1033 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1034 value += valuesBuffer[i];
1035 }
1036
1037 value *= G4UniformRand();
1038
1039 i = n;
1040
1041 while (i > 0)
1042 {
1043 --i;
1044
1045 if (valuesBuffer[i] > value)
1046 {
1047 delete[] valuesBuffer;
1048 return i;
1049 }
1050 value -= valuesBuffer[i];
1051 }
1052
1053 if (valuesBuffer)
1054 delete[] valuesBuffer;
1055
1056 }
1057 } else
1058 {
1059 G4Exception("G4DNARuddIonisationModel::RandomSelect",
1060 "em0002",
1062 "Model not applicable to particle type.");
1063 }
1064
1065 return level;
1066}
1067
1068//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1069
1070G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1071{
1072 G4double sigma = 0.;
1073
1074 const G4DynamicParticle* particle = track.GetDynamicParticle();
1075 G4double k = particle->GetKineticEnergy();
1076
1077 G4double lowLim = 0;
1078 G4double highLim = 0;
1079
1080 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1081
1082 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1083 pos1 = lowEnergyLimit.find(particleName);
1084
1085 if (pos1 != lowEnergyLimit.end())
1086 {
1087 lowLim = pos1->second;
1088 }
1089
1090 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1091 pos2 = highEnergyLimit.find(particleName);
1092
1093 if (pos2 != highEnergyLimit.end())
1094 {
1095 highLim = pos2->second;
1096 }
1097
1098 if (k >= lowLim && k <= highLim)
1099 {
1100 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1101 pos = tableData.find(particleName);
1102
1103 if (pos != tableData.end())
1104 {
1105 G4DNACrossSectionDataSet* table = pos->second;
1106 if (table != 0)
1107 {
1108 sigma = table->FindValue(k);
1109 }
1110 } else
1111 {
1112 G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1113 "em0002",
1115 "Model not applicable to particle type.");
1116 }
1117 }
1118
1119 return sigma;
1120}
1121
1122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1123
1124G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1125 const G4String& /* particle */)
1126{
1127 return 0;
1128}
1129
G4AtomicShellEnumerator
@ eIonizedMolecule
G4double S(G4double temp)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:180
G4double G4Log(G4double x)
Definition: G4Log.hh:227
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
const double C2
#define C1
#define G4UniformRand()
Definition: Randomize.hh:52
static G4Alpha * Alpha()
Definition: G4Alpha.cc:88
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)
static G4DNAGenericIonsManager * Instance(void)
G4ParticleDefinition * GetIon(const G4String &name)
const std::vector< G4double > * GetNumMolPerVolTableFor(const G4Material *) const
Retrieve a table of molecular densities (number of molecules per unit volume) in the G4 unit system f...
static G4DNAMolecularMaterial * Instance()
G4ParticleChangeForGamma * fParticleChangeForGamma
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *p, G4double ekin, G4double emin, G4double emax)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
G4DNARuddIonisationModel(const G4ParticleDefinition *p=0, const G4String &nam="DNARuddIonisationModel")
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
Definition: G4Electron.cc:93
static G4LossTableManager * Instance()
G4VAtomDeexcitation * AtomDeexcitation()
const G4Material * GetMaterial() const
size_t GetIndex() const
Definition: G4Material.hh:255
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:691
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
G4double GetPDGCharge() const
const G4String & GetParticleName() const
Definition: G4Pow.hh:49
static G4Pow * GetInstance()
Definition: G4Pow.cc:41
G4double logZ(G4int Z) const
Definition: G4Pow.hh:137
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:162
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:230
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:87
static G4Proton * Proton()
Definition: G4Proton.cc:92
const G4DynamicParticle * GetDynamicParticle() 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:746
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:600
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:124
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:802
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:607
void ProposeTrackStatus(G4TrackStatus status)
const G4Track * GetCurrentTrack() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
const G4double pi