Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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
39//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
40
41using namespace std;
42
43//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44
46 const G4String& nam) :
47G4VEmModel(nam), isInitialised(false)
48{
49 fpWaterDensity = 0;
50
51 slaterEffectiveCharge[0] = 0.;
52 slaterEffectiveCharge[1] = 0.;
53 slaterEffectiveCharge[2] = 0.;
54 sCoefficient[0] = 0.;
55 sCoefficient[1] = 0.;
56 sCoefficient[2] = 0.;
57
58 lowEnergyLimitForZ1 = 0 * eV;
59 lowEnergyLimitForZ2 = 0 * eV;
60 lowEnergyLimitOfModelForZ1 = 100 * eV;
61 lowEnergyLimitOfModelForZ2 = 1 * keV;
62 killBelowEnergyForZ1 = lowEnergyLimitOfModelForZ1;
63 killBelowEnergyForZ2 = lowEnergyLimitOfModelForZ2;
64
65 verboseLevel = 0;
66 // Verbosity scale:
67 // 0 = nothing
68 // 1 = warning for energy non-conservation
69 // 2 = details of energy budget
70 // 3 = calculation of cross sections, file openings, sampling of atoms
71 // 4 = entering in methods
72
73 if (verboseLevel > 0)
74 {
75 G4cout << "Rudd ionisation model is constructed " << G4endl;
76 }
77
78 // Define default angular generator
80
81 // Mark this model as "applicable" for atomic deexcitation
83 fAtomDeexcitation = 0;
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;
132 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
133 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
134 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
135 G4ParticleDefinition* 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 G4DNAGenericIonsManager *instance;
303
304 if (
305 particleDefinition != G4Proton::ProtonDefinition()
306 &&
307 particleDefinition != instance->GetIon("hydrogen")
308 &&
309 particleDefinition != instance->GetIon("alpha++")
310 &&
311 particleDefinition != instance->GetIon("alpha+")
312 &&
313 particleDefinition != instance->GetIon("helium")
314 )
315
316 return 0;
317
318 G4double lowLim = 0;
319
320 if ( particleDefinition == G4Proton::ProtonDefinition()
321 || particleDefinition == instance->GetIon("hydrogen")
322 )
323
324 lowLim = lowEnergyLimitOfModelForZ1;
325
326 if ( particleDefinition == instance->GetIon("alpha++")
327 || particleDefinition == instance->GetIon("alpha+")
328 || particleDefinition == instance->GetIon("helium")
329 )
330
331 lowLim = lowEnergyLimitOfModelForZ2;
332
333 G4double highLim = 0;
334 G4double sigma=0;
335
336 G4double waterDensity = (*fpWaterDensity)[material->GetIndex()];
337
338 const G4String& particleName = particleDefinition->GetParticleName();
339
340 // SI - the following is useless since lowLim is already defined
341 /*
342 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
343 pos1 = lowEnergyLimit.find(particleName);
344
345 if (pos1 != lowEnergyLimit.end())
346 {
347 lowLim = pos1->second;
348 }
349 */
350
351 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
352 pos2 = highEnergyLimit.find(particleName);
353
354 if (pos2 != highEnergyLimit.end())
355 {
356 highLim = pos2->second;
357 }
358
359 if (k <= highLim)
360 {
361 //SI : XS must not be zero otherwise sampling of secondaries method ignored
362
363 if (k < lowLim) k = lowLim;
364
365 //
366
367 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
368 pos = tableData.find(particleName);
369
370 if (pos != tableData.end())
371 {
372 G4DNACrossSectionDataSet* table = pos->second;
373 if (table != 0)
374 {
375 sigma = table->FindValue(k);
376 }
377 }
378 else
379 {
380 G4Exception("G4DNARuddIonisationModel::CrossSectionPerVolume","em0002",
381 FatalException,"Model not applicable to particle type.");
382 }
383
384 }
385
386 if (verboseLevel > 2)
387 {
388 G4cout << "__________________________________" << G4endl;
389 G4cout << "G4DNARuddIonisationModel - XS INFO START" << G4endl;
390 G4cout << "Kinetic energy(eV)=" << k/eV << " particle : " << particleDefinition->GetParticleName() << G4endl;
391 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
392 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
393 //G4cout << " - Cross section per water molecule (cm^-1)="
394 //<< sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
395 G4cout << "G4DNARuddIonisationModel - XS INFO END" << G4endl;
396 }
397
398 return sigma*waterDensity;
399
400}
401
402//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
403
404void G4DNARuddIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
405 const G4MaterialCutsCouple* couple,
406 const G4DynamicParticle* particle,
407 G4double,
408 G4double)
409{
410 if (verboseLevel > 3)
411 {
412 G4cout << "Calling SampleSecondaries() of G4DNARuddIonisationModel"
413 << G4endl;
414 }
415
416 G4double lowLim = 0;
417 G4double highLim = 0;
418
419 G4DNAGenericIonsManager *instance;
421
422 if ( particle->GetDefinition() == G4Proton::ProtonDefinition()
423 || particle->GetDefinition() == instance->GetIon("hydrogen")
424 )
425
426 lowLim = killBelowEnergyForZ1;
427
428 if ( particle->GetDefinition() == instance->GetIon("alpha++")
429 || particle->GetDefinition() == instance->GetIon("alpha+")
430 || particle->GetDefinition() == instance->GetIon("helium")
431 )
432
433 lowLim = killBelowEnergyForZ2;
434
435 G4double k = particle->GetKineticEnergy();
436
437 const G4String& particleName = particle->GetDefinition()->GetParticleName();
438
439 // SI - the following is useless since lowLim is already defined
440 /*
441 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
442 pos1 = lowEnergyLimit.find(particleName);
443
444 if (pos1 != lowEnergyLimit.end())
445 {
446 lowLim = pos1->second;
447 }
448 */
449
450 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
451 pos2 = highEnergyLimit.find(particleName);
452
453 if (pos2 != highEnergyLimit.end())
454 {
455 highLim = pos2->second;
456 }
457
458 if (k >= lowLim && k <= highLim)
459 {
460 G4ParticleDefinition* definition = particle->GetDefinition();
461 G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
462 /*
463 G4double particleMass = definition->GetPDGMass();
464 G4double totalEnergy = k + particleMass;
465 G4double pSquare = k*(totalEnergy+particleMass);
466 G4double totalMomentum = std::sqrt(pSquare);
467 */
468
469 G4int ionizationShell = RandomSelect(k,particleName);
470
471 G4double bindingEnergy = 0;
472 bindingEnergy = waterStructure.IonisationEnergy(ionizationShell);
473
474 //SI: additional protection if tcs interpolation method is modified
475 if (k<bindingEnergy) return;
476 //
477
478 // SI - For atom. deexc. tagging - 23/05/2017
479 G4int Z = 8;
480 //
481
482 G4double secondaryKinetic = RandomizeEjectedElectronEnergy(definition,k,ionizationShell);
483
484 G4ThreeVector deltaDirection =
485 GetAngularDistribution()->SampleDirectionForShell(particle, secondaryKinetic,
486 Z, ionizationShell,
487 couple->GetMaterial());
488
489 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),deltaDirection,secondaryKinetic);
490 fvect->push_back(dp);
491
492 // Ignored for ions on electrons
493 /*
494 G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
495
496 G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
497 G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
498 G4double finalPz = totalMomentum*primaryDirection.z() - deltaTotalMomentum*deltaDirection.z();
499 G4double finalMomentum = std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
500 finalPx /= finalMomentum;
501 finalPy /= finalMomentum;
502 finalPz /= finalMomentum;
503
504 G4ThreeVector direction;
505 direction.set(finalPx,finalPy,finalPz);
506
507 fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
508 */
509
511
512 // sample deexcitation
513 // here we assume that H_{2}O electronic levels are the same of Oxigen.
514 // this can be considered true with a rough 10% error in energy on K-shell,
515
516 size_t secNumberInit = 0;// need to know at a certain point the energy of secondaries
517 size_t secNumberFinal = 0;// So I'll make the diference and then sum the energies
518
519 G4double scatteredEnergy = k-bindingEnergy-secondaryKinetic;
520
521 // SI: only atomic deexcitation from K shell is considered
522 if(fAtomDeexcitation && ionizationShell == 4)
523 {
524 const G4AtomicShell* shell
525 = fAtomDeexcitation->GetAtomicShell(Z, G4AtomicShellEnumerator(0));
526 secNumberInit = fvect->size();
527 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, 0, 0);
528 secNumberFinal = fvect->size();
529
530 if(secNumberFinal > secNumberInit)
531 {
532 for (size_t i=secNumberInit; i<secNumberFinal; ++i)
533 {
534 //Check if there is enough residual energy
535 if (bindingEnergy >= ((*fvect)[i])->GetKineticEnergy())
536 {
537 //Ok, this is a valid secondary: keep it
538 bindingEnergy -= ((*fvect)[i])->GetKineticEnergy();
539 }
540 else
541 {
542 //Invalid secondary: not enough energy to create it!
543 //Keep its energy in the local deposit
544 delete (*fvect)[i];
545 (*fvect)[i]=0;
546 }
547 }
548 }
549
550 }
551
552 //This should never happen
553 if(bindingEnergy < 0.0)
554 G4Exception("G4DNAEmfietzoglouIonisatioModel1::SampleSecondaries()",
555 "em2050",FatalException,"Negative local energy deposit");
556
557 //bindingEnergy has been decreased
558 //by the amount of energy taken away by deexc. products
559 if (!statCode)
560 {
563 }
564 else
565 {
568 }
569
570 // debug
571 // k-scatteredEnergy-secondaryKinetic-deexSecEnergy = k-(k-bindingEnergy-secondaryKinetic)-secondaryKinetic-deexSecEnergy =
572 // = k-k+bindingEnergy+secondaryKinetic-secondaryKinetic-deexSecEnergy=
573 // = bindingEnergy-deexSecEnergy
574 // SO deexSecEnergy=0 => LocalEnergyDeposit = bindingEnergy
575
576 // TEST //////////////////////////
577 // if (secondaryKinetic<0) abort();
578 // if (scatteredEnergy<0) abort();
579 // if (k-scatteredEnergy-secondaryKinetic-deexSecEnergy<0) abort();
580 // if (k-scatteredEnergy<0) abort();
581 /////////////////////////////////
582
583 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
585 ionizationShell,
586 theIncomingTrack);
587 }
588
589 // SI - not useful since low energy of model is 0 eV
590
591 if (k < lowLim)
592 {
596 }
597
598}
599
600//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
601
602G4double G4DNARuddIonisationModel::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition,
603 G4double k,
604 G4int shell)
605{
606 G4double maximumKineticEnergyTransfer = 0.;
607
608 G4DNAGenericIonsManager *instance;
610
611 if (particleDefinition == G4Proton::ProtonDefinition()
612 || particleDefinition == instance->GetIon("hydrogen"))
613 {
614 maximumKineticEnergyTransfer = 4. * (electron_mass_c2 / proton_mass_c2) * k;
615 }
616
617 else if (particleDefinition == instance->GetIon("helium")
618 || particleDefinition == instance->GetIon("alpha+")
619 || particleDefinition == instance->GetIon("alpha++"))
620 {
621 maximumKineticEnergyTransfer = 4. * (0.511 / 3728) * k;
622 }
623
624 G4double crossSectionMaximum = 0.;
625
626 for (G4double value = waterStructure.IonisationEnergy(shell);
627 value <= 5. * waterStructure.IonisationEnergy(shell) && k >= value;
628 value += 0.1 * eV)
629 {
630 G4double differentialCrossSection =
631 DifferentialCrossSection(particleDefinition, k, value, shell);
632 if (differentialCrossSection >= crossSectionMaximum)
633 crossSectionMaximum = differentialCrossSection;
634 }
635
636 G4double secElecKinetic = 0.;
637
638 do
639 {
640 secElecKinetic = G4UniformRand()* maximumKineticEnergyTransfer;
641 }while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition,
642 k,
643 secElecKinetic+waterStructure.IonisationEnergy(shell),
644 shell));
645
646 return (secElecKinetic);
647}
648
649//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
650
651// The following section is not used anymore but is kept for memory
652// GetAngularDistribution()->SampleDirectionForShell is used instead
653
654/*
655 void G4DNARuddIonisationModel::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition,
656 G4double k,
657 G4double secKinetic,
658 G4double & cosTheta,
659 G4double & phi )
660 {
661 G4DNAGenericIonsManager *instance;
662 instance = G4DNAGenericIonsManager::Instance();
663
664 G4double maxSecKinetic = 0.;
665
666 if (particleDefinition == G4Proton::ProtonDefinition()
667 || particleDefinition == instance->GetIon("hydrogen"))
668 {
669 maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
670 }
671
672 else if (particleDefinition == instance->GetIon("helium")
673 || particleDefinition == instance->GetIon("alpha+")
674 || particleDefinition == instance->GetIon("alpha++"))
675 {
676 maxSecKinetic = 4.* (0.511 / 3728) * k;
677 }
678
679 phi = twopi * G4UniformRand();
680
681 // cosTheta = std::sqrt(secKinetic / maxSecKinetic);
682
683 // Restriction below 100 eV from Emfietzoglou (2000)
684
685 if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
686 else cosTheta = (2.*G4UniformRand())-1.;
687
688 }
689 */
690//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
691G4double G4DNARuddIonisationModel::DifferentialCrossSection(G4ParticleDefinition* particleDefinition,
692 G4double k,
693 G4double energyTransfer,
694 G4int ionizationLevelIndex)
695{
696 // Shells ids are 0 1 2 3 4 (4 is k shell)
697 // !!Attention, "energyTransfer" here is the energy transfered to the electron which means
698 // that the secondary kinetic energy is w = energyTransfer - bindingEnergy
699 //
700 // ds S F1(nu) + w * F2(nu)
701 // ---- = G(k) * ---- -------------------------------------------
702 // dw Bj (1+w)^3 * [1 + exp{alpha * (w - wc) / nu}]
703 //
704 // w is the secondary electron kinetic Energy in eV
705 //
706 // All the other parameters can be found in Rudd's Papers
707 //
708 // M.Eugene Rudd, 1988, User-Friendly model for the energy distribution of
709 // electrons from protons or electron collisions. Nucl. Tracks Rad. Meas.Vol 16 N0 2/3 pp 219-218
710 //
711
712 const G4int j = ionizationLevelIndex;
713
714 G4double A1;
715 G4double B1;
716 G4double C1;
717 G4double D1;
718 G4double E1;
719 G4double A2;
720 G4double B2;
721 G4double C2;
722 G4double D2;
723 G4double alphaConst;
724
725 // const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
726 // The following values are provided by M. dingfelder (priv. comm)
727 const G4double Bj[5] = { 12.60 * eV, 14.70 * eV, 18.40 * eV, 32.20 * eV, 540
728 * eV };
729
730 if (j == 4)
731 {
732 //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
733 A1 = 1.25;
734 B1 = 0.5;
735 C1 = 1.00;
736 D1 = 1.00;
737 E1 = 3.00;
738 A2 = 1.10;
739 B2 = 1.30;
740 C2 = 1.00;
741 D2 = 0.00;
742 alphaConst = 0.66;
743 } else
744 {
745 //Data For Liquid Water from Dingfelder (Protons in Water)
746 A1 = 1.02;
747 B1 = 82.0;
748 C1 = 0.45;
749 D1 = -0.80;
750 E1 = 0.38;
751 A2 = 1.07;
752 // Value provided by M. Dingfelder (priv. comm)
753 B2 = 11.6;
754 //
755 C2 = 0.60;
756 D2 = 0.04;
757 alphaConst = 0.64;
758 }
759
760 const G4double n = 2.;
761 const G4double Gj[5] = { 0.99, 1.11, 1.11, 0.52, 1. };
762
763 G4DNAGenericIonsManager* instance;
765
766 G4double wBig = (energyTransfer
767 - waterStructure.IonisationEnergy(ionizationLevelIndex));
768 if (wBig < 0)
769 return 0.;
770
771 G4double w = wBig / Bj[ionizationLevelIndex];
772 // Note that the following (j==4) cases are provided by M. Dingfelder (priv. comm)
773 if (j == 4)
774 w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
775
776 G4double Ry = 13.6 * eV;
777
778 G4double tau = 0.;
779
780 G4bool isProtonOrHydrogen = false;
781 G4bool isHelium = false;
782
783 if (particleDefinition == G4Proton::ProtonDefinition()
784 || particleDefinition == instance->GetIon("hydrogen"))
785 {
786 isProtonOrHydrogen = true;
787 tau = (electron_mass_c2 / proton_mass_c2) * k;
788 }
789
790 else if (particleDefinition == instance->GetIon("helium")
791 || particleDefinition == instance->GetIon("alpha+")
792 || particleDefinition == instance->GetIon("alpha++"))
793 {
794 isHelium = true;
795 tau = (0.511 / 3728.) * k;
796 }
797
798 G4double S = 4. * pi * Bohr_radius * Bohr_radius * n
799 * std::pow((Ry / Bj[ionizationLevelIndex]), 2);
800 if (j == 4)
801 S = 4. * pi * Bohr_radius * Bohr_radius * n
802 * std::pow((Ry / waterStructure.IonisationEnergy(ionizationLevelIndex)),
803 2);
804
805 G4double v2 = tau / Bj[ionizationLevelIndex];
806 if (j == 4)
807 v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
808
809 G4double v = std::sqrt(v2);
810 G4double wc = 4. * v2 - 2. * v - (Ry / (4. * Bj[ionizationLevelIndex]));
811 if (j == 4)
812 wc = 4. * v2 - 2. * v
813 - (Ry / (4. * waterStructure.IonisationEnergy(ionizationLevelIndex)));
814
815 G4double L1 = (C1 * std::pow(v, (D1))) / (1. + E1 * std::pow(v, (D1 + 4.)));
816 G4double L2 = C2 * std::pow(v, (D2));
817 G4double H1 = (A1 * std::log(1. + v2)) / (v2 + (B1 / v2));
818 G4double H2 = (A2 / v2) + (B2 / (v2 * v2));
819
820 G4double F1 = L1 + H1;
821 G4double F2 = (L2 * H2) / (L2 + H2);
822
823 G4double sigma =
824 CorrectionFactor(particleDefinition, k) * Gj[j]
825 * (S / Bj[ionizationLevelIndex])
826 * ((F1 + w * F2)
827 / (std::pow((1. + w), 3)
828 * (1. + G4Exp(alphaConst * (w - wc) / v))));
829
830 if (j == 4)
831 sigma = CorrectionFactor(particleDefinition, k) * Gj[j]
832 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
833 * ((F1 + w * F2)
834 / (std::pow((1. + w), 3)
835 * (1. + G4Exp(alphaConst * (w - wc) / v))));
836
837 if ((particleDefinition == instance->GetIon("hydrogen"))
838 && (ionizationLevelIndex == 4))
839 {
840 // sigma = Gj[j] * (S/Bj[ionizationLevelIndex])
841 sigma = Gj[j] * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
842 * ((F1 + w * F2)
843 / (std::pow((1. + w), 3)
844 * (1. + G4Exp(alphaConst * (w - wc) / v))));
845 }
846
847 // if ( particleDefinition == G4Proton::ProtonDefinition()
848 // || particleDefinition == instance->GetIon("hydrogen")
849 // )
850
851 if (isProtonOrHydrogen)
852 {
853 return (sigma);
854 }
855
856 if (particleDefinition == instance->GetIon("alpha++"))
857 {
858 slaterEffectiveCharge[0] = 0.;
859 slaterEffectiveCharge[1] = 0.;
860 slaterEffectiveCharge[2] = 0.;
861 sCoefficient[0] = 0.;
862 sCoefficient[1] = 0.;
863 sCoefficient[2] = 0.;
864 }
865
866 else if (particleDefinition == instance->GetIon("alpha+"))
867 {
868 slaterEffectiveCharge[0] = 2.0;
869 // The following values are provided by M. Dingfelder (priv. comm)
870 slaterEffectiveCharge[1] = 2.0;
871 slaterEffectiveCharge[2] = 2.0;
872 //
873 sCoefficient[0] = 0.7;
874 sCoefficient[1] = 0.15;
875 sCoefficient[2] = 0.15;
876 }
877
878 else if (particleDefinition == instance->GetIon("helium"))
879 {
880 slaterEffectiveCharge[0] = 1.7;
881 slaterEffectiveCharge[1] = 1.15;
882 slaterEffectiveCharge[2] = 1.15;
883 sCoefficient[0] = 0.5;
884 sCoefficient[1] = 0.25;
885 sCoefficient[2] = 0.25;
886 }
887
888 // if ( particleDefinition == instance->GetIon("helium")
889 // || particleDefinition == instance->GetIon("alpha+")
890 // || particleDefinition == instance->GetIon("alpha++")
891 // )
892 if (isHelium)
893 {
894 sigma = Gj[j] * (S / Bj[ionizationLevelIndex])
895 * ((F1 + w * F2)
896 / (std::pow((1. + w), 3)
897 * (1. + G4Exp(alphaConst * (w - wc) / v))));
898
899 if (j == 4)
900 sigma = Gj[j]
901 * (S / waterStructure.IonisationEnergy(ionizationLevelIndex))
902 * ((F1 + w * F2)
903 / (std::pow((1. + w), 3)
904 * (1. + G4Exp(alphaConst * (w - wc) / v))));
905
906 G4double zEff = particleDefinition->GetPDGCharge() / eplus
907 + particleDefinition->GetLeptonNumber();
908
909 zEff -= (sCoefficient[0]
910 * S_1s(k, energyTransfer, slaterEffectiveCharge[0], 1.)
911 + sCoefficient[1]
912 * S_2s(k, energyTransfer, slaterEffectiveCharge[1], 2.)
913 + sCoefficient[2]
914 * S_2p(k, energyTransfer, slaterEffectiveCharge[2], 2.));
915
916 return zEff * zEff * sigma;
917 }
918
919 return 0;
920}
921
922//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
923
924G4double G4DNARuddIonisationModel::S_1s(G4double t,
925 G4double energyTransferred,
926 G4double slaterEffectiveChg,
927 G4double shellNumber)
928{
929 // 1 - e^(-2r) * ( 1 + 2 r + 2 r^2)
930 // Dingfelder, in Chattanooga 2005 proceedings, formula (7)
931
932 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
933 G4double value = 1. - G4Exp(-2 * r) * ((2. * r + 2.) * r + 1.);
934
935 return value;
936}
937
938//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
939
940G4double G4DNARuddIonisationModel::S_2s(G4double t,
941 G4double energyTransferred,
942 G4double slaterEffectiveChg,
943 G4double shellNumber)
944{
945 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 2 r^4)
946 // Dingfelder, in Chattanooga 2005 proceedings, formula (8)
947
948 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
949 G4double value = 1.
950 - G4Exp(-2 * r) * (((2. * r * r + 2.) * r + 2.) * r + 1.);
951
952 return value;
953
954}
955
956//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
957
958G4double G4DNARuddIonisationModel::S_2p(G4double t,
959 G4double energyTransferred,
960 G4double slaterEffectiveChg,
961 G4double shellNumber)
962{
963 // 1 - e^(-2 r) * ( 1 + 2 r + 2 r^2 + 4/3 r^3 + 2/3 r^4)
964 // Dingfelder, in Chattanooga 2005 proceedings, formula (9)
965
966 G4double r = R(t, energyTransferred, slaterEffectiveChg, shellNumber);
967 G4double value = 1.
968 - G4Exp(-2 * r)
969 * ((((2. / 3. * r + 4. / 3.) * r + 2.) * r + 2.) * r + 1.);
970
971 return value;
972}
973
974//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
975
976G4double G4DNARuddIonisationModel::R(G4double t,
977 G4double energyTransferred,
978 G4double slaterEffectiveChg,
979 G4double shellNumber)
980{
981 // tElectron = m_electron / m_alpha * t
982 // Dingfelder, in Chattanooga 2005 proceedings, p 4
983
984 G4double tElectron = 0.511 / 3728. * t;
985 // The following values are provided by M. Dingfelder (priv. comm)
986 G4double H = 2. * 13.60569172 * eV;
987 G4double value = std::sqrt(2. * tElectron / H) / (energyTransferred / H)
988 * (slaterEffectiveChg / shellNumber);
989
990 return value;
991}
992
993//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
994
995G4double G4DNARuddIonisationModel::CorrectionFactor(G4ParticleDefinition* particleDefinition,
996 G4double k)
997{
998 G4DNAGenericIonsManager *instance;
1000
1001 if (particleDefinition == G4Proton::Proton())
1002 {
1003 return (1.);
1004 } else if (particleDefinition == instance->GetIon("hydrogen"))
1005 {
1006 G4double value = (std::log10(k / eV) - 4.2) / 0.5;
1007 // The following values are provided by M. Dingfelder (priv. comm)
1008 return ((0.6 / (1 + G4Exp(value))) + 0.9);
1009 } else
1010 {
1011 return (1.);
1012 }
1013}
1014
1015//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1016
1017G4int G4DNARuddIonisationModel::RandomSelect(G4double k,
1018 const G4String& particle)
1019{
1020
1021 // BEGIN PART 1/2 OF ELECTRON CORRECTION
1022
1023 // add ONE or TWO electron-water ionisation for alpha+ and helium
1024
1025 G4int level = 0;
1026
1027 // Retrieve data table corresponding to the current particle type
1028
1029 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1030 pos = tableData.find(particle);
1031
1032 if (pos != tableData.end())
1033 {
1034 G4DNACrossSectionDataSet* table = pos->second;
1035
1036 if (table != 0)
1037 {
1038 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
1039
1040 const size_t n(table->NumberOfComponents());
1041 size_t i(n);
1042 G4double value = 0.;
1043
1044 while (i > 0)
1045 {
1046 i--;
1047 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
1048 value += valuesBuffer[i];
1049 }
1050
1051 value *= G4UniformRand();
1052
1053 i = n;
1054
1055 while (i > 0)
1056 {
1057 i--;
1058
1059 if (valuesBuffer[i] > value)
1060 {
1061 delete[] valuesBuffer;
1062 return i;
1063 }
1064 value -= valuesBuffer[i];
1065 }
1066
1067 if (valuesBuffer)
1068 delete[] valuesBuffer;
1069
1070 }
1071 } else
1072 {
1073 G4Exception("G4DNARuddIonisationModel::RandomSelect",
1074 "em0002",
1076 "Model not applicable to particle type.");
1077 }
1078
1079 return level;
1080}
1081
1082//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1083
1084G4double G4DNARuddIonisationModel::PartialCrossSection(const G4Track& track)
1085{
1086 G4double sigma = 0.;
1087
1088 const G4DynamicParticle* particle = track.GetDynamicParticle();
1089 G4double k = particle->GetKineticEnergy();
1090
1091 G4double lowLim = 0;
1092 G4double highLim = 0;
1093
1094 const G4String& particleName = particle->GetDefinition()->GetParticleName();
1095
1096 std::map<G4String, G4double, std::less<G4String> >::iterator pos1;
1097 pos1 = lowEnergyLimit.find(particleName);
1098
1099 if (pos1 != lowEnergyLimit.end())
1100 {
1101 lowLim = pos1->second;
1102 }
1103
1104 std::map<G4String, G4double, std::less<G4String> >::iterator pos2;
1105 pos2 = highEnergyLimit.find(particleName);
1106
1107 if (pos2 != highEnergyLimit.end())
1108 {
1109 highLim = pos2->second;
1110 }
1111
1112 if (k >= lowLim && k <= highLim)
1113 {
1114 std::map<G4String, G4DNACrossSectionDataSet*, std::less<G4String> >::iterator pos;
1115 pos = tableData.find(particleName);
1116
1117 if (pos != tableData.end())
1118 {
1119 G4DNACrossSectionDataSet* table = pos->second;
1120 if (table != 0)
1121 {
1122 sigma = table->FindValue(k);
1123 }
1124 } else
1125 {
1126 G4Exception("G4DNARuddIonisationModel::PartialCrossSection",
1127 "em0002",
1129 "Model not applicable to particle type.");
1130 }
1131 }
1132
1133 return sigma;
1134}
1135
1136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1137
1138G4double G4DNARuddIonisationModel::Sum(G4double /* energy */,
1139 const G4String& /* particle */)
1140{
1141 return 0;
1142}
1143
G4AtomicShellEnumerator
@ eIonizedMolecule
double S(double temp)
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:179
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
const double C2
#define C1
#define G4UniformRand()
Definition: Randomize.hh:52
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: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)
G4double GetPDGCharge() const
const G4String & GetParticleName() const
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: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 ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
const G4double pi